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ABSTRACT 

The BOXPLX of a user specified cost function is evaluated 
as a control system design method. The emphasis is on compen- 
sator design for multiple-input-multiple-output plants, but the 
technique should be applicable to any control system with 
adjustable parameters. 

The engineer must define the plant dynamics, inputs and 
outputs, the desired output responses for the specified inputs, 
weighting factors for the performance measure, the type of com- 
pensation to be used, the free parameters, and initial values 
and constraints on the free parameters. 

The program gives a direct search of the optimum values, 
within the specified constraints, for the free parameters. 
Evaluation of the results is given in graphical plots of the 
output time responses including the desired output response as 
well as the compensated response. 

The Transfer Matrix Method is used to provide an analytical 
check on the accuracy of the method and the procedure is illus- 
trated with a two-input-two-output system example. 



4 



TABLE OF CONTENTS 



I . INTRODUCTION 8 

II. MULTIVARIABLE SYSTEM COMPENSATION 10 

A. TRANSFER MATRIX-ANALYTICAL METHODS 10 

1. Introduction 10 

2. Feedback Compensation 11 

3. Cascade Compensation 14 

4. Cascade and Feedback Compensation 16 

3. PARAMETER OPTIMIZATION - COMPUTER METHOD 19 

1. Optimization Technique 19 

2. Performance Measure 20 

C. A TWO-INPUT-TWO-OUTPUT SYSTEM EXAMPLE 26 

1. Transfer Matrix Solution 26 

2. Parameter Optimization Solution 32 

III. CONCLUSIONS 5 2 

APPENDIX A Computer Program 53 

LIST OF REFERENCES 7 8 

INITIAL DISTRIBUTION LIST 79 



LIST OF FIGURES 



2-1 Multivariable Plant 10 

2-2 Generalized Multivariable Plant 10 

2-3 Multivariable Plant With Feedback Compensation 12 

2-4 Matrix/Vector Representation of Multivariable Plant 

With Feedback Compensation 14 

2-5 Multivariable Plant with Cascade Compensation 15 

2-6 Matrix/Vector Representation of Multivariable Plant 

With Cascade Compensation 16 

2- 7 Matrix/Vector Representation of Multivariable Plant 

With Cascade and Feedback Compensation 16 

3- 1 Multivariable Plant to be Compensated 27 

3-2 Multivariable Plant with Generalized Compensation 27 

3-3 Compensated Multivariable Plant 33 

3-4A Desired Output/Actual Output Yl & Z1 vs Time 

(Uncompensated) R1 = Unit Step; R2 = 0 3 4 

3-4B Yl and Z1 vs Time with Rl = Unit Step, R2 = 0.0 

For Compensated System 3 4 

3-5A Y2 and Z2 vs Time with Rl = Unit Step, R2 = 0.0 

(Y2 = -3.57 X lO"^) 35 

ss 

3-5B Y2 and Z2 vs Time with Rl = Unit Step, R2 = 0.0 

For Compensated System 35 

3-6A Yl and Zl vs Time with Rl = 0.0, R2 = Unit Step 

(Yl = 2.91 X 10“^) 36 

ss 

3-6B Yl and Zl vs Time with Rl = 0.0, R2 = Unit Step 

For Compensated System 3 6 

3-7A Y2 and Z2 vs Time With Rl = 0.0, R2 = Unit Step 37 

3-7B Y2 and Z2 vs Time With Rl = 0.0, R2 = Unit Step 

For Compensated System 3 7 

3-8A Yl and Zl vs Time with Rl = Unit Step, R2 = Unit 

Step 3 8 



6 



38 



3-8B Yl and Z1 vs Time with R1 = Unit Step, R2 = Unit 

Step For Compensated System 

3-9A Y2 and Z2 vs Time with Rl = Unit Step, R2 = Unit 

Step 39 

3-9B Y2 and Z2 vs Time with Rl = Unit Step, R2 = Unit 

Step for Compensated System 39 

3-lOA Yl and Reference 1 vs Time Simulated with Thesis 

Program Parameters Set at Analytical Optimiom 42 

3-lOB Y2 and Reference 2 vs Time Simulated with Thesis 

Program Parameters -Set at' Analytical Optimum 43 

3-llA Yl and Reference 1 vs Time After Optimization 

With Integral-Square Error Performance Measure 44 

3-llB Y2 and Reference 2 vs Time After Optimization 

With Integral Square-Error Performance i!4easure 
(Curve With Greatest Overshoot is ZZ) 45 

3-12A Yl and Reference 1 vs Time After Optimization 

With Time Multiplied Integral-Square Error 
Performance Measure 46 

3-12B Y2 and Reference 2 vs Time After Optimization 

With Time Multiplied Integral Square-Error 
Performance Measure 47 

3-13A Yl and Reference 1 vs Time After Optimization 

With Integral Absolute Error Performance 
Measure 43 

3-13B Y2 and Reference 2 vs Time After Optimization 

With Integral Absolute Error Performance 
Measure 4 9 

3-14A Yl and Reference 1 vs Time After Optimization 

With Time Multiplied Integral Absolute Error 
Performance Measure 50 

3-14B Y2 and Reference 2 vs Time After Optimization 

With Time Multiplied Integral Absoluts Error 
Performance Measure 



51 



I. 



INTRODUCTION 



As control systems become more and more complex, design 
engineers are faced with an increasingly difficult problem in 
the development of compensators. Single-input-single-output 
linear systems can be compensated using Bode plots or root 
locus plots and trial and error methods to meet the system 
specifications [1] . When compensating multiple-input-multiple- 
output linear systems, use can be made of transfer matrix 
techniques described by Ogata [2] to achieve total decoupling 
and the desired output responses. Since decoupling during 
the transient period of real systems may require controls 
which exceed physical constraints imposed by the system, it is 
sometimes desirable to decouple during steady state conditions 
and allow coupling during the transient period [ 3] . 

Non-linear, single-input, single-ouput systems, wherein the 
non-linearity can be lumped into one element, can be analyzed 
and compensated effectively using the phase plane method (for 
second order systems) and describing function analysis [2,4] . 

Some of these nvethods require extensive mathematical manipulations 
and become rather unwieldy, but do produce the desired results. 

In all of the analysis and design approaches mentioned above, 
any one of several well established simulation programs can be 
used to evaluate the design. 

To reduce the tedium involved in the design of compensators, 
frequency domain and time domain computer optimization schemes 
have been developed, which will set free parameters in the 
compensator to yield a system response which most nearly 
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approximates a reference response specified by the user. Lima 
used this approach to address the problem of controller design 
for ships engaged in underway replenishment [5 ] . MacNamara 
applied this method to an autopilot design Csl , and Vines 
developed a generalized program for compensator optimization 
in single-input-single-output systems [7]. 

The following chapters present, a generalized discussion of 

/ 

transfer matrix analysis and compensation of multivariable 
systems, optimization techniques, and performance measures. 

A two- input- two-output system is then compensated using transfe 
matrix techniques, and using parameter optimization. The 

/ 

results are then compared. 
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II. MULTIVARIABLE SYSTEM COMPENSATION 



A. TRANSFER J^ATRIXtANALYTICAL METHODS 
1. Introduction 

A linear, multiple-input-multiple-output or multivariable 
system can be described by a transfer matrix, which is simply 
an extension of the transfer function concept. Consider the 
two-input-two-output open loop plant: 



R1(S) 

R2(Sl 



Gp(S) 



CKS) 



Figure 2-1 Multivariable Plant 



The Block Diagram of the plant can be represented in the 
general sense as shown in Figure 2-2. 




Figure 2-2 Generalized Multivariable Plant 
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Algebraic equations for the outputs in the independent variable 
S, can be written directly: 



Cl(s) = GplKs)Rl(s) + Gpl2(s) R2{s) 
C2(s) = Gp21(s)Rl{s) + Gp22(s) R2(s) 



( 1 ) 

( 2 ) 



If the inputs Rj (s) and the outputs C^(s) 
as vectors of one dimension (i, j = 1, 2) 
written in matrix form: 



Cl(s) 




Gpll (s) 


Gpl2 (s) 


_C2(s)_ 




_Gp21(s) 


Gp22 (s)_ 



are each considered 
the equations can be 



RKs) 

R2{s) 



( 3 ) 



The number of inputs and outputs can be increased, and need not 
be equal so that in the general case: 



Cl(s) 
C2 (s) 
1 



1 

1 

Cj_(s) 



or : 







~ — -r 


Gpll(s) Gpl2(s) Gpij (s) 




■Rl(s) 


Gp21 (s) 




R2(s) 


1 






1' 






1 






Gpil (s) Gpij (s) 




R. (s) 


— , — 




J 



C (s) = Gp (s) R(s) 



(4) 



(5) 



where Gp(s) is the transfer matrix (ixj) of the plant. 



2. Feedback Compensation 

Transfer matrices can be used to describe the compen- 
sation of multivariable systems as well as the plant. Consider 
the two-input-two-ouput plant above with feedback compensation 
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as shown in Figure 2-3. 




Figure 2-3 

Multivariable Plant with Feedback Compensation 



Algebraic equations for R^(s) and R2 (s) can be written: 



R1(S) = U1(S) 


- H ll(s) 


CKs) 


- Hl2(s) 


C2(s) 


(6) 


R2(s) = U2(s) 


- H 21(s) 


CKs) 


- H22(s) 


C2 (s) 


(7) 



with C, R, and U two element, 
equations can be rewritten: 



— 










Rl(s) 




Ul(s) 




HlKs) 


R2 (s) 




U2 (s) 




H21 (s) 



one dimension vectors these 







— 1 


H12 (s) 




CKs) 


H22 (sj^ 




_C2(s)_ 



12 



or in the general case: 



P.(s) = U(s) - H(s)C(s) (9) 

where H(s) is the transfer matrix (jxi) of the compensation. 
From the uncompensated plant: 

C(s) = G (s)R(s) (10) 

ir 

Substituting equation (9) for R(s) in Equation (10) : 

C(s) = G (s) [ U(s) - H(s) C(s)l (11) 

~P ...... 

Rearranging : 

-[I + G (s)H(s)]C(s) = G (s)U(s) (12) 

.. -P - -p - 

Since Gp(s) is an ixj matrix, K(s) is a jxi matrix the product, 
Gp(s)H(s), is a square, ixi matrix. Assigning that [I + Gp(s) 
is non-singular the inverse can be taken. 

Premultiplying both sides of equation (12) by [I+G^ ( s) H ( s ) ] 

C(s) = [I + G^ (s)H(s) (s)U(s) (13) 

- .. ..p .. ,p , 

C(s) = G(s)U(s) (14) 

The compensated system transfer matri.x (ixj) is then: 

G(s) = [I + G^(s)H(s) ] ’^G„(S) (15) 

The plant with feedback compensation can be depicted more 
clearly as shown in Figure 2-4. 
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I m 



R[S) 




Figure 2-4 

Matrix/Vector Representation of Multivariable 
Plant with Feedback Compensation 

To determine the feedback compensation, H(s), required to 

achieve design specifications, the dynamics of the plant must 

be known and put into transfer matrix form, G (s) . The desired 

~? 

closed loop dynamics must then be reflected in the closed loop 
transfer matrix, G(s) . Then through matrix manipulation the 
compensation, H(s), can be computed in closed form provided 
that G(s) and G^(s) are square and non-singular. 

Premultiplying both sides of equation (15) by [I+G (s)H(s)] 
and then postmultiply ing by G(s)~^: 

I + G^(s)H(s) = G^(s)G(s)"^ (16) 

- ~p - ..p 

Rearranging : 

G (s)H(s) = G (s)G(s)"^ - I (17) 

p ^ ^ 

Then multiplying both sides of equation (16) by G (s)~^ 

H(s) = G(s)”^ - Gp(s)"^ (13) 

3. Cascade Compensation 

Cascade compensation can be treated in a similar fashion. 
Consider the same basic plant with cascade compensation as shown 



14 



in Figure 2-5. 




Figure 2-5 

.Multivariable Plant with Cascade Compensation 

The algebraic equations for R^(s) and are: 

Rl(s) = Gcll(s)U^(s) + Gcl2(s)U2(s) 

R2(s) = Gc21(s)Ul(s) + Gc22(s)U2(s) 



Rl(s) 




"gcII (s) 


Gc 1 2 ( s ) 




"ui (s) 


Jl2(s)_ 




_Gc21 (s) 


Gc22 (s)_ 




JJ2 (s) _ 



in the general case: 

C(s) = Gp (s) Gc (s) U(s) 

G(s) = Gp(s)Gc(s) 

If Gp(s) is a square non-singular matri:< the compensation 
Gc(s) is given by: 

Gc(s) = Gp(s)"^ G(s) 



(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 
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The compensated plant can be depicted as shown in Figure 2-6, 



RlSl 




U(S) 




C(Sl 






Figure 2-6 

Matrix/Vector Representation of Multivariable 
Plant with Cascade Compensation 



4. Cascade and Feedback Compensation. 

Cascade and feedback compensation can also be combined 
as shown in Figure 2-7. 




Figure 2-7 

Matrix/Vector Representation of Multivariable 
Plant with Cascade and Feedback Compensation 
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E(s) = 


R(s) - H(s)C(s) 


(25) 




U(s) = 


G^ (s) E (s) 


(26) 


C(s) 


= G^(s)U(s) 


= Gp(s)Gc(s) [R(s) - H(s)C(s)] 


(27) 


C(s) 


= Gp(s)Gc(s)R(s) - Gp(s)Gc(s)H(s)C(s) 


(28) 


Gp (s) Gc (s) R(s) = 


[I + Gp(s)Gc(s)H(s)]C(s) 


(29) 



If Gp(s) is an ixj matrix Gc(s) will be a j x j matrix, and 
H(s) will be jxi, in which case Gp (s) Gc (s) H (s) is a square 

matrix (ixi) . If [I + 'Go (s) Gc ( s) His ) ] is non-singular, it can be 

inverted and: 

C(s) = [I + Gp(s) Gc (s) H (s) ] “^Gp (s) Gc (s)R (s) (30) 

C(s) — G(s) R(s) 

The compensated system transfer matrix G(s) (ixj) is : 

G(s) = [I + Gp (s) GC (s) H (s) ] Gp(s)Gc(s) (31) 

Equation (31) contains two unknown matrices, H(s) and Gc(s) . 

One of these must be selected arbitrarily and then in certain 



special cases 


the other can 


be solved explicitly. 


For example 




H(s) = 


[; ;] ■ < 


(32) 


In which case: 








G ( s ) = 


[I + Gp(s)Gc(s)] ^ Gp(s)Gc(s) 


(33) 


[I + Gp(s)Gc(s) ]G(s) = 


Gp ( s ) Gc ( s ) 


(34) 


G(s) = 


Gp(s)Gc(s) [I 


- G ( s ) ] 


(35) 


If Gp(s) is a 


square matrix 


G(s) and Gc(s) will be 


square. If 


Gp(s) [G(s)"^ 


- I] , and G(s) 


are non-singular; 
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Gp (s) 


-1 


= Gc(s) [G(s) - I] 


(36) 


Gc (s) 




= Gp(s)“^ - 1]“^ 


(37) 


If Gc(s) 


and 


{[G(s) ^ - l]Gp(s)} are assximed to be 


non-singular 


equation 


(37) 


can be rearranged for convenience : 




Gp (s) 


Gc (s 


) = [G(s)"^ - I]"^ 


(38) 


Gp(s) 


= 


[G(s)"^ - 1]“^ Gc(s)'"^ 


(39) 


Gc(s) 


-1 


= [G(s)“^ - I] Gp(s) 


(40) 


Gc (s) 


= 


{ [G(s)“^ - I] Gp(s) 


(41) 



This result can also be obtained directly from equation (37) 
using the matrix identity (AB) ^ A~^. 

In each of the systems discussed above the compensation 
cannot be solved for in closed form, except in special cases 



as indicated. 
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B. PARAMETER OPTIMIZATION - COMPUTER METHOD 



1. Optimization Technique 

The Complex Method of constrained optimization (BoxPLX) 
was utilized in this program- The algorithm suggested by 
Box (9) was modified by the programmer, R. R. Hilleary, Naval 
Postgraduate School, to find the constrained minimum, and 
includes an integer programming option, as suggested in (10) . 

The Complex Method is a modification of the unconstrained, direct 
search. Simplex Method introduced in (11). Direct search methods 
compare function and constraint values only in searching for 
the minimum value of the function. These direct search methods 
have been widely used because of their robustness, reliability, 
and ease in programming and use. They have widespread applic- 
ability. The one drawback is that they are generally less 
efficient than gradiant-based techniques (12). A more detailed 
discussion of the Complex Method is included in the program 
documentation. 

In this program, the main program simulates the reference 
output and then calls FUNCTION BOXPLX which in turn calls FUNCTION 
FE. FUNCTION FE calls SUBROUTINE Plant which simulates the 
compensated system. FE then computes the performance measure 
and returns to BOXPLX. BOXPLX keeps track of the compensator 
parameter values, computes a new set of free parameters and 
calls FE again. This iteration continues until termination as 
described in the documentation. See Appendix .A for more 
information on the program. 
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2. Performance Measures 



In order for a numerical method for parameter opti- 
mization to produce a result, a decision process must be com- 
pleted, the result of which is the set of parameters which is 
defined as ''optimum.” This term optimum is an enigima of sorts. 
The program addressed here utilizes a comparison of time domain 
response of a compensated system with a reference response 
specified by the user. The word optimum above is dependent on 
the designer's knowledge of what reference response is best 
suited to the output he is trying to control. Furthermore, 
with the exception of very simple linear systems, it is diffi- 
cult to match a reference response exactly, so the designer 
must decide if the comparison at steady state is more or less 
important than that during the transient period. Additionally 
very high frequencies with periods less than one half the 
integration step size may be filtered out by the digital process 
in the computer, and not be present in the simulated output. 

The designer is faced with some decisions which he must make 
objectively before the ”optim;m'' set of parameters can be 
computed . 

The decision process mentioned above is based upon a 
performance measure which is an indicator of the "goodness" 
of a given set of parameters being varied in the optimization 
process. The designer must select the performance measure; the 
performance measure is then minimized by numerical mechods , 
since the objective here is to have the system response match 
the reference response. 
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The term selectivity is used to describe a quality of the 
performance measures. If three dimensional space is considered 
with the performance measure plotted as a function of two free 
variables in a compensator of the system, the resultant topology 
should include one or more minima for the performance measure. 
The minima with the lowest value are referred to as the global 
minima, although the minimization performed is constrained. 
Selectivity is directly related to the steepness of the slope 
of the contour in the area immediately adjacent to the minima. 
This selectivity is a function of the performance measure used 
as well as the plant and compensator dynamics (2,3) . 

Several performance measures have been suggested in numerous 
sources, some of which will be discussed here. In each case 
the performance measure will be represented by the variable J. 

Consider a single- input single-output system with output, 
c(t) , and reference response, r(t) . The error, e(t) , for this 
case is defined as the difference between the reference res- 
ponse and the system response or: 

e(t) = r(t) - c(t) (42) 



Integral square-error. The defining equation for finite 



time is : 



J = f e^(t)dt 

This performance measure is not very selective and tends 
emphasize large errors and ignore small ones (2,8). 



(43) 



to 
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Integral of time multiplied square-error . The defining 
equation for finite time is: 



■r 



t e (t)dt 



(44) 



Steady state error is penalized more heavily by this per- 
formance measure. This performance measure is said to be 
more selective than the integral square-error criterion in that 
the value of the integral tends to change more rapidly with 
changes in the system parameters (2,8). 

Integral absolute error . The defining equation for finite 
time is : 

rt 

J = le(t)ldt (45) 

Jo 

This performance measure is slightly more selective than the 
integral square-error criterion (2,8). 

Integral of time-multiplied absolute error . The defining 
equation for finite time is: 



J 




1 1 e ( t) 1 dt 



(46) 



This performance measure is more selective than those mentioned 
previously and, as in the case with the integral of time multi- 
plied square-error criterion, this performance measure tends to 
emphasize errors late in the transient period and i.n steady 
state and deemphasize those which occur early in the transient 
period (2,3) . 
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The optimization of the performance measure can be addressed 
in two categories; the case where the reference response speci- 
fied is the forcing function; or the reference response is some 
function, usually a second order response, which the designer 
believes to be a model response for his application. In the 
first case the peak overshoot, setting time, rise time, and 
natural frequency, are determined by the dynamics of the com- 
pensated system and by the performance measure selected. How- 
ever, if the designer desires a second order response with a 
certain minimal rise time, he can specify a second order 
reference response which meets this requirement, and use the 
integral absolute error or integral error-squared performance 
criteria to optimize the compensation parameters. 

Provisions for both of these approaches have been included 
in the program in that the user can reference a second order 
response, where he specifies the damping factor, natural frequency 
and gain, or he can utilize the table loop-up feature in which 
he can specify 'any reference response he desires. 

Other performance measures can be developed using classical 
measures of system performance such as maximum overshoot, time 
for the error to reach its first zero, time to reach maximum 
overshoot, settling time, frequency of oscillation of the 
transient or steady state error. However, performance measures 
based on these characteristics would be complex and increase 
computation time, and their desirability over those previously 
discussed is questionable. 

The performance measure used in the optimization program 
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(see Appendix) is the sum of the weighted performance measures 
of each output of the multivariable system. 



M 



J 




(47) 



The subscript i denotes the output for which J. ^ is calculated. 
The subscript k denotes one of four performance measures which 
can be selected by the user to be applied to each output in 
turn : 

Integral square error 

N 



Integral of time multiplied square error 

N 

Ji,2 = S ( (Rj ,i - C. , ) 
j=l 3.1 

Integral of absolute-error 

N 



Ji,l 



Z 

j = x 



(Rj,i - C ) 
J > 



2 



(48) 



J i , 3 = Z I R j , i 

j=l 




(50) 



Integral o 



i. 



time multiplied square-error 



Ji , 4 



N 

S I Rj ,i - C. J (iAT) 
J t ^ 

j = l 



(51) 



N represents the number of integration steps in the si.mulation , 
AT represents the integration step size, .M is the number of 
outputs being optimized, and W is the weighting factor applied 



1 





to each output. 

It should be noted that the performance measures are an 
approximation and when the weighting factor is introduced 
the performance measure diverges significantly from the defining 
integral equation. Comparison of performance measures for a 
given plant with various types of compensation should be done 
with the same weighting factors. The same performance measure 
could be achieved for a system using two different compensation 
schemes, since the performance measures for each channel may 
vary significantly. 
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C. A TWO INPUT-'I^'JO OUTPUT SYSTEM EXAMPLE 

1. Transfer Matrix Solution 

In order to demonstrate the program it was decided to 
analyze a simple system, synthesize the compensation required 
to produce a specified output, and then show that the optimi- 
zation program will arrive at the same results. It is emphasized 
that this program will not determine the type of compensation 
required to achieve the desired response. It will set free 
parameters in the compensators (specified by the user) to most 
nearly produce the desired response. 

The analysis and compensation was performed using the 
transfer matrix method described by OGATA (2). The uncompensated 
plant and the compensated plant were simulated to show that ohe 
compensation achieved the design specifications- Certain para- 
meters in the compensator were defined as free parameters. 

These free parameters were offset from their opti.mal values 
and constraints were introduced which restricted the range of 
values for these parameters during optimization. The plant 
and compensators with free parameters were simulated in the 
optimization program and optimization was accomplished. 

To illustrate, a simple two-input two-output linear 
system shown in Figure 3-1 was selected as the plant to be 
compensated. The design specifications chosen were: 

1. Channel one and channel two are to be decoupled during 
the transient period as well as in steady state. 

2. Channel one is to have a second order response to 

a step input with a natural frequency, of 10, a 
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Figure 3-1 Multivariable Plant to be Compensated 




Figure 3-2 Multivariable Plant with Generalized Compensation 
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m 






damping factor, of 0.4, and no steady state error. 

3. Channel two is to have a second order response to a 
step input with a natural frequency, W^, of 4, a damping 
factor, C / of 0.6, and no steady state error. 

To achieve total decoupling the off-diagonal elements of 
the transfer matrix of the compensated system must be zero. 
The desired compensated system transfer matrix, G(s), is: 



G(s) 



100 

+ 8S + 100 

0 



16 



S + 4.8S + 16 



(52) 



The input/output equations obtained from Figure 3-1 are: 



Y^(s) 



U^(s) 



(S+7) (S + 12; 



Y2(s) 



(s+12) 



(53) 



Y2(3) 



^2^^^ (s+2) (s+4.3) 



Y^(s) 



4 

(s + 2) 



(54) 



Substitution of equation (53) into equation (54), and equation 
(54) into equation (53) produces equations (55) and (56) . 

U^(s) 2U2 (s) 3Y^(s) 

^1^^^ ^ (s+7) (s + 12) ^ (s + 2) (s+4.3) (s + 12) ^ (s + 2) (s + 12) 

U2(s) 4Uj^(s) 8Y2 (s) 

^2^^^ "" (s+2) (s+4.3) ^ (s+12) (s+7) (s+2) ^ (s+2) (s+12) 



Collecting like terms in equations (55) and (56) and rearranging 
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2U^ (s) 



( 57 ) 



„ s^+14s+16 _ ^1^^^ 

(s+2) (s+12) (s+7) (s+12) (s+2) (s+4.3) (s+12) 



4U3_(s) 



Y (s) s +14s+16 

^2^^^ (s+2) (s+12) (s+12) (s+7) (s+2) 



U2(s) 



(s+2) (s+4.3) 



(58) 



Multiplying both sides of equations (57) and (53) by 

(s+2) (s+12) 
s^+14s+16 



Y^(s) = 



(s+2) (s) 



2U^ (s) 



(s+7) (s +14S+16) (s+4.3) (s +14s+16) 



(59) 



Y2(s) = 



4U^(s) 



(s+12)U2(s) 



2 2 

(s+7) (s +14S+16) (s+4.3) (s +14s+16) 



(60) 



The transfer Matrix for the uncompensated plant is: 



G^(s) 




(61) 



(s+7) (s^+14s+16) (s+4.3) (2^+14s+16) 



The generalized compensation to be used is shown in Figure 
3-2. Equation (41) will be used to calculate the required 
compensation; since in this case the feedback ccmpensarion 
transfer matrix, H(s), is the identity matrix. 



G (s) = {[G(s)’^ - I]G (s)}"-^ 



(41) 



Since G (s) is a 2x2 matrix in this example G (s) and G(s) will 
also be 2x2. 3y inspection (Equation (52)) , G(s) is non-singular 
and can be inverted. 
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The compensated system transfer matrix G(s) is given in 
equation (52). The inverse of G(s) is calculated by equation 
(62) . 



G(s) 



-1 



COF 

G(s) 



(62) 



Where COF is the transpose of the cofactor matrix of G(s) , 
and |G(s)| is the determinant of G(s). 

16 0 



COF‘ 



s +4.8s+16 



100 



s +8S+100 



(63) 



G(s) = 



(16) (IGO) 



(s^+3s+100) (s^+4.8s+16) 



(64) 



VJ ^ a ) - 



s i-8s+100 
100 



’ + 4 . 3 s+1 6 

Te 



(65) 



'“s^+8s 
100 



[G(s)~^-I] = 



s ■(•4.8s 
16 



( 66 ) 



The plant transfer matrix, G^(s) , is given in equation (61) as 
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G ( s ) — 

- 



(s+2) 2 

(s+7) (s^+14s+16) (s+4.3) (s^+14s+16) 



(s+12) 



(s+7) (s^+14s+16) (s+4.3) (s^+14s+16) 



(61) 



Multiplying equation (61) by equation (66) : 



[ G(s)"^-I]G^(s) 

- - - ~ p ~ 



s (s+8) (s+2) 



2s (s+8 ) 



100 (s+7) (s^+14s+16) 100(s+4.3) (s^+14s+16) 



(67) 



4 (s+4 .8) 



s(s+4.8) (s+12) 



16(s+7) (s^+14s+16) 16(s+4.3) (s^+14s+16) 



If this matrix is non-singular it can be inverted 



_ 1 



[G ( s ) —I] G^ ( 3 ) 



s (s+8) (S+2) (s+4. 8) (s + 12) 

fl6) flOO) (s+7) (s+4.3) (s^+14s+16)^ 



( 68 ) 



3 b (;34-4.8) (s+8) 



(16) (100) (s+7) (s+4.3) (s^+14sfl6)^ 



1 [G(s)“^ - I 



s^(s+4) (s+8) (s^+14s + 24-8) 

(16) (100) (s+7) (s+4.3) (s^+14s+16)^ 

s^(s+4,8) (s+8) 

(16) (100) (s+7) (s+4.3) (s^ + 14s + 16) 



(69) 
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T 

COF 



sCs + 4.8) (S-H2) 
16(s+4.3) (s^+14s+16) 

-4s (s+4.8) 

16 (s+7) (s^+14s+16) 



-2s(s + 8) 

100(s + 4.3) (s^ + 14s + 16) 

s (s+8) (s+2) 

100 (s+7) (s^+14s+16) 



(70) 



G^(s) = { [G(s)“^ - l]G„(s)} 



100 (s+7) (s+12) 
s (s+8) 



-32 (s+7) 
s (s+4 . 8 ) 



-400(s+4.3) 
s ( s+8) 



16 (s+2) (s+4. 3) 
s(s+4.3) 



(71) 



One of the possible realizations of this compensation is 
shown in Figure 3-3. To obtain a confirmation of the design, 
the system was simulated with and without compensation using 
International Business Machines' Digital Simulation Language, 

DSL (13). Figures 3-4 through 3-9 show the results of that 
simulation, which indicate that the compensators used do in 
fact achieve the stated objectives of the design. In each plot, 
trace number 1 is the desired response, Z, and trace number 2 
the system response, Y. In Figures 3-4B, 3-73, 3-83, and 3-9B 
the two traces are superimposed. In Figures 3-5A and 3-6A 
the desired responses are zero. The reader should note the 
scale factors present on some of these plots. 



2. Parameter Optimization Solution 

The compensated system and the desired response curves 



were then simulated and solutions were obtained using the 
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Ficjure 3-3 Compensated Multivaiiable Plant 
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Y1 and Zl vs Time with Rl = Unit Step, 




Y1 and Zl vs Time with Rl = Unit Step, 
For Compensated System 
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R2 = 0.0 



t 



R2 = 0.0 
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Figure 3-5A 

Y2 and Z2 vs Time with Rl = Unit Step, R2 = 0.0 

(Y2 = 3.57 X 10“^) 
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Figure 3-5B 

Y2 and Z2 vs Time with Rl = Unit Step, R2 = 0.0 
For Compensated System 
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Figure 3-6A 

Y1 and Z1 vs Time with Rl = 0.0, R2 = Unit Step 

(Yl^^ = 2.91 X lO"^) 
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Figure 3-6B 

Y1 and Z1 vs Time with Rl = 0.0, R2 = Unit Step 

For Compensated System 
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Y2 and Z2 vs Time with R1 = 0.0, R2 = Unit Step 

For Compensated System 
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Figure 3-8A 

Y1 and Zl vs Time with Rl = Unit Step, R2 = Unit Step 
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Figure 3-8B 

Y1 and Zl vs Time with Rl = Unit Step, R2 

For Compensated System 



Unit Step 
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Figure 3-9A 

Y2 and Z2 vs Time with Rl = Unit Step, R2 = Unit Step 




Y2 and Z2 vs Time with Rl = Unit Step, R2 - Unit Step 

For Compensated System 
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program discussed in this section. The time response results 
are shown in Figure 3-10, There is a slight difference between 
the desired response curve and the output of the simulated 
system. 

Next, two parameter optimization was performed. The pro- 
portional gain constant of each compensator was selected as 
the adjustable parameters. Each of the four performance measures 
given by equations (47) through (51) were used in successive 
optimizations . All other aspects of the program were unchanged 
during this optimization. Equal weighting was applied to each 
channel. The gains 100 and 16 were chosen because their 
variation should have similar impact on the two respective out- 
puts. The upper limits were 105 and 21 and the lower limits 
were 95 and 11 respectively. The integration step size was 
.02 and the final time was 5 seconds. Optimization was begun 
with the two gains set at 105 and 11 respectively, and Rl = R2 = 
Unit Step functions. The results are given in Table 3-1. 

The system was simulated using these values for the two 
free parameters. The simulation results are shown in Figures 
3-11 through 3-14. This simulation indicates the time multiplied 
integral square-error performance measure produced the best: 
results , 
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Performance 

Measure 



Value of 

Minimum *^100 

Performance 

Measure 



G 



16 



Integral , 

Square- 1.2026097 x 10"'^ 104.9999 20.99998 

Error 



Time 

Multiplied c 

Integral 1.641641 x 10“ 102.2554 15.27091 

Square- 

Error 



Integral - 

Absolute 1.353197 x lo"-^ 104.9862 14.82275 

Error 



Time 

Multiplied . 

Integral 9.994698 x 10~ 103.2728 14.81198 

Absolute 

Error 



Table 3-1 



RESULTS OF OPTIMIZATION 
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Figure 3-lOA 

Y1 and Reference 1 vs Time Simulated with Thesis Program 
Parameters Set at Analytical Optimum 
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Figure 3-lOB 

Y2 and Reference 2 vs Tine Simulated with Thesis Program 
Parameters Set at Analytical Optimum 
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Figure 3-llA 

Yl and Reference 1 vs Tine After Optinization With 
Integral-Square Error Performance Measure 
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Figure 3-113 

Y2 and Reference 2 vs Time After Optimization With 
Integral Square-Error Performance Measure 
(Curve With Greatest Overshoot is ZZ) 
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Figure 3-12A 

Yl and Reference 1 vs Time After Optimization With 
Time Multiplied Integral-Square Error Performance Measure 
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Figure 3-12B 

Y2 and Reference 2 vs Time After Optimization With 
Time 2iultiplied Integral Square-Error Performance Measure 
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Figure 3-13A 

Yl and Reference 1 vs Time After Optimization With 
Integral Absolute Error Performance Measure 
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Figure 3-13B 

Y2 and Reference 2 vs Time After Optimization With 
Integral Absolute Error Performance Measure 



49 



XIO 

out) 002 Oou 006 006 OtO 012 



Y1 

Z1 



3 * 

a 




t 

4i___ 

'Ctoa 002 004 006 CQ8 010 

Figure 3-14A 

Y1 and Reference 1 vs Time After Optimization With 
Time Multiplied Integral Absolute Error Performance Measure 
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Figure 3-143 

Y2 and Reference 2 vs Time After Optimization With 
Time Multiplisd Integral Absolute Error Psrformance Measure 
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CONCLUSIONS 



From the results of the two-input two-output example 
system, the following conclusions can be stated: 

1. The technique of cost function minimization in the 
time domain can be used effectively for compensator parameter 
optimization in multiple input multiple output control systems. 

2. Cost function minimization is not a sxibstitute for 
control system design. The user must be able to select the 
proper type of compensator which is capable of achieving the 
desired results, before initiating the optimization of the 
free parameters . 

3. The four performance measures evaluated produced 
slightly different solutions to the same problem, with the 
time multiplied integral square-error yielding the best 
results in the example given. 
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APPENDIX A 



COMPENSATOR OPTIMIZATION IN 
MULTIVARIABLE CONTROL SYSTEMS 
BY 



JOHN T. MOWREY 
MAJOR USMC 
DECEMBER 1978 

THIS PROGRAM WAS DEVELOPED AS A GENERAL PURPOSE DESIGN 
AID TO BE USED tq OPTIMIZE FREE =>ARAMETERS IN THE COMP- 
ENSATORS OF A MULTIVARIABLE CONTROL SYSTEM, IT CAN ALSO 
8 6 USED FOR. THE OPTIMIZATION OF' FREE PARAMETERS IN A 
SINGLE INPUT SINGLE OUTPUT SYSTEM' OP ^0 SIMULATE EITHER 
TYPE OF SYSTEM. OPTIMIZATION IS ACCOMPLISHED BY MINIMIZ- 
ING THE ERROR BETWEEN A USER SPECIFIED REFERENCE CURVE 
AND THE OUTPUT OF THE SIMULATED COMPENSATED SYSTEM. 



DATA CAROS 

:«t T)t ****:*** # 

CARO 1 NRUNStNaPT,NGRAPri,NREF,IIN,NOLT,IPRINT,IFREC 
FORMAT! 815) 

NRUNS-FLAG FOR NUMBER OF RUNS. NRUNS MUST EQUAL 
1 WHEN OPTIMIZING. 

NOPT-FLAG FOR OPT IM IZAT ION / S I TULA TI ON 

1- OPTIMIZAT ION 

2- SIMULATION ONLY 

(IF N0PT=1 THE plant WILL BE SIMULATED AFTER 
optimization AND THE REQUESTED OUTPUT PRODUCED 

NGRAPH-FLAG FOR TY°E OF GRAPHICAL OUTPUT D5SIPEC 

1- PRODUCES PRINTER plQtS 

2- oRO DUCES VERSA TEC PLOTS 

3 - NO GRAPHS PRODUCED 

NREF-fLAG FOR REFERENCE CURVE 

1 - PRO DUCES CURVE 

2- NO REFERENCE CURVE IS PRODUCED 

IIN-THE number of inputs (UP tq 25) 

NOUT-THE NUMBER OF OUTPUTS IU° TO 3) 

I PRINT-FLAG FOR TABULATED DATA FOR REFERENCE 
CURVE AND SYSTEM RESPONSE. 

1- PRQDUCES PRINTED OUTPUT 

2- NO OUTPUT DATA PRINTED 

IPREQ-CREOUENCY OF 3RINTED OUTPUT. IF IFREQ=10 
6VF?Y lOTH DATA POINT /^ILL BE PRINTED. 



CARD 2 T47,DT,TF 

FORMAT(3riO,5 ) 
TA7-INITIAL TIME 
DT-INT6GRATI0N STEP SIZE 
TF-FINAL time 
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CARO 3 NV ,NAV ,NTA,,NPRf IP, I PM 
F3RMAT( 615 ) 

THIS CARD MUST 3E OMITTED IF N0PT=2 

NV-THE NUMBER OF FREE VARIABLES BE SET BY 
OPTIMIZATION. 

NAV-THE NUMBER OF AUXILLIARY VARIABLES 

NTA-THE NUMBER OF TRIALS ALLOWED BOXPLX 

NPR-FREQUENCY OF OUTPUT FROM BOXPLX FOR 
DIAGNOSTIC PURPOSES .. 

IP-FLAG FOR INTEG£R/REAL*3 OPTIMIZATION. 

^ * 

^ NOTE * 

* * 

SEE 80XPLX FOR ADDITIONAL INFORMATION ON THE 5 VARIABLES 
LISTED ABOVE. 30XPLX WAS CONVERTED TO DOUBLE PRECISION 
FOR USE IN THIS PROGRAM; HOV»£VER THE DOCUMENTATION WAS 
NOT CHANGED 

IPM-CLAG FOR PERFORMANCE MEASURE. 

I'HE performance MEASURES USED ARE PROPORT lONAL 
TO THOSE INDICATED BELOW. 

1- INTEGRAL SQUARE-ERROR (WEIGHTED) 

2- TIME MULTIPLIED INTEGRAL SQUARE-ERROR 
(WE IGHT ED ) 

3- IN^EGRAL ABSOLUTE ERROR (WEIGHTED) 

4- TIME MULTIPLIED PlTEGRAL ABSOLUTE ERROR 
(WEIGHTED) 



CAROS 4-7 XS(I ) 

FORM-AT(8FiO. 5) 

ALL OF these CAROS mijst BE QmiTTEO I= M0PT=2 
XS(I) ARE THE STARTING VALUES OF THE FREE PARA- 
METERS. THE NUMBER OF CARDS REQUIRED IS DEPEN- 
DENT ON NV. ONE VALLE MUST BE READ IN FCR EACH 
FREE PARAMETER. IE IF NV= B ONLY CARO 4 IS RE- 
QUIRED AND 3 PARAMETER VALUES SHOULD BE ON IT 
IN 8F10, 5 FORMAT. 



CARDS 8-11 XU(I) 

format (8F10.5) 

THESE CARDS ARE THE SAME AS CAROS 4-7 EXCEPT 
XU (I) ARE THE UPPER LIMITS OM THE FREE PARA- 
METERS. DMI T IF N0PT=2 . 



CARDS 12-15 XL( I) 

FORM AT (3 FiO .5 ) 

these cards are the same as CAROS 4-7 EXCEPT 
XL(I) ARE THE LOWER LIMITS ON THE FREE PARA- 
METERS. OMIT IF NCPT =2. 
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c 

c 

n 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



CAROS 16»18, AND 20 I OUT » I WT , I TA3 

=ORMAT( 3F10.5) 

I OUT-THE 'NUMBSR OF THE BLOCK FROM WHICH THE OUT- 
PUT IS TAKEN, I F 2 OR MORE BLOCKS FEED AN OUT- 
PUT NODE, A TYOE 1 BLOCK WITH G=1 MUST BE 
PLACED BETWEEN THAT NODE ANC THE OUTPUT. 

FOR A 3 OUTPUT SYSTEM ALL 3 CARDS ARE REQUIRED, 

FOR A 2 OUTPUT SYSTEM 15 AND 13 ARE REQUIRED. 

IWT-AN INTEGER WEIGHTING FUNCTION WHICH WILL 
85 CONVERTED TO REAL=^8.AMD WHICH ALLOWS THE 
USER TO penalize THE ERROR BETWEEN THE SYSTEM 
OUTPUT AND THE REFERENCE RESOONSE MORE HEAVILY 
AT ONE OUTPUT THAN ANOTHER. INTEGERS BETWEEN 1 
AND 10 ARE RECOMMENDED. WEIGHTING SHOULD BE 
REDUCED IF OVERFLOW IS ENCOUNTERED. 

Itab-A flag for TABLE LOOK UP OF THE REFERENCE 
C urves. 

1- OROGRAM will READ TABULATED DATA FOR THE 
REFERENCE CURVES. 

2- PROGRAM WILL CALCULATE A SECOND ORDER RES- 
PONSE CURVE. 



CAROS 17t19, and 21-B ET A , DELTA , WN , AM© ,DEL AY , I NP UT 
=ORMAT( 5F10.5 ,AA) 



C 

r 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



these cards input the data required to COMOU’E a 

SECOND ORDER RESPONSE CURVE. 1= THE USER DESIRES 
the table look up feature these CAROS ARE RE- 
PLACED BY TABULATED DATA. IF A TABULATED REF- 
ERENCE CURVE IS DESIRED FOR OLTPUT NC. 1, CARD 
17 IS REPLACED BY TABULATED DATA IN cFlO. 5 FORMAT 
the number of data points REQUIRED IS T F/ OT (TF 
AND DT FROM CARD 2). 



cm 



B ETA ( I ) 



R ( I ) 



( S**2 + 2^0ELTA( I ) *WN( I )*S^■WN ( I )*’f'2) 



cm IS THE DESIRED REFERENCE RESPONSE ASSOC- 
IATED WITH OUTPUT lOUT(I). 

Rm IS THE FORCING FUNCTION S = ECIFIEC BY: 
AM0(I)-THE amplitude of the FORCING FUNCTICN 
OELAYd )-0ELAY AFTER T 47 BEFORE THE FORCING 
FUNCTION IS APOLIEDITIVE COMAIN). 
INPUT(I)-TH5 TYPE OF FORCING F UNCTION : STE=> , 
RAMP, OR PARA, 

3ETA(I)-THE gain of the TPaNSFcR FUNCTION. 
OELTAm-THE DAMPING FACTOR. 

WN (I ) -THE NATURAL FREQUENCY. 

CARD 22 N 

FORMAT! I2J 

THE NUMBER OF BLOCKS IN THE SYSTEM (UP TO 25). 
CARDS 23-47 IC ( I ) , ID ( I ) , IE ( I ) , I V ( I ) , G < I ) , p ( I ) , Z { I ) 



FORMAT! 4I2t 2X,3F1C.5) 

THE NUMBER OF CAROS REQUIRED IS DETER'^INED BY N. 
IC!I»-THE NUMBER OF THE BLOCK! 1-25). 
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ID(I)-THE TYPE 3F BL3CK. SIX TYPES ARE AVAItaBLE 



ID ( I ) =1 



1 E ( I ) 



IV ( I) 



lOU ) =2 



I£( I ) 



S+P 1 



IV( I ) 



lom =3 



IE ( I ) 



IV{ I ) 



I 0 ( I ) =4 




IE(I ) I 



ion ) =5 




IV(I)=G*IE{n FOR G*IE(I) LE U. 

AMD G^I£( I ) GE LL 

IV( I) =UL FOR G*IE(I ) GT LL 

IV( I ) = LL FOR G*IE( II LT UL 



IO ( I ) =6 



I£( I) 



/ -8 
// j G 



/I 

B / I 



IV( I) 



IV(I)=G^I£(I) =OR G*I£(I) GE S 
OR LE -8 

IV(I)=0 FOR G»I£{I) LT 3 AND GT -3 
IE(I)-THE INPUT NODE NUMBER. 

IV(I)-THE OUTPUT NODE NUMBER. 

G(I)-TH£ GAIN FIELD. G(I) REPRESENTS THE GAIN 
OF THE 3L0CK. 

P( I )-THE POLE FIELD. 

FOR ID( I ) =2 OP 3 P(I ) = c>OlE 
F-^R IO(Il=<+ ?{ n=OAMPING FACTOR 
FOR IO(I)=5 P{n = UP'RER LIMIT, UL. 

=OR I0(I)=6 P( I ) =REFERENCE, 8. 

U I )-th£ zero field 

FOR 10(1 ) =3 Z( I ) =ZERO 

FOR I0(I)=4 Z(I» = NATUPAL FO EOU ENC Y , WN . 

FOR I0(I)=5 Z(1)=LCWER LImIT,lL. 
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c 

CAROS 48-72 AMP ( n , 0ELAY( I ) , IMPUT( n , lEEd ) C 

F0RMAT(2F10.5tA4, 12) C 

C 

ONE CARO IS NEEDED FOR EACH FORCING FLMCTION ON C 
THE SIMULATED COMPENSATED SYSTEM. C 

C 

AMP( n, 0ELAY( I )» AN0INPUT( I ) ARE AS DESCRIBED FOR C 
CAROS 17, 19, AND 21, EXCEPT THEY APPLY TC THE COM- C 
PENSATEO SYSTEM IN THIS CASE. C 

C 

IEE(I)-THE NODE TC )»HICH THE FORCING FUNCTION IS C 
APPLIED. C 

C 

CAROS 73-78 TITLE CAROS FOR VERSATEC FLCTS. IF NGRAPH=2, C 
NOUT*2 TITLE CARDS ARE REQUIRED(2 PER OUTPUT) C 

TITLES GO IN COLUMNS 1-48 ONLY, C 

" C 

* * C 

* NOTE ^ C 

3S * C 

WHEN CPTI Mj ZING , ONE FORTRAN CARD ®ER PREF VARIABLE MUST C 



3E INSERTED IN SUBROUTINE PLANT. THE APORQPRIATE POSITIQNC 
IS IDENTIFIED BY THE COMMENT CARO ' ENTER G(I)=C(I) VAL- C 
UE S A T TH I S POI NT • . C 

C 

C 

iik-sie :^3ae :iK ie lie ^ :ir :ie :be lieia: lieiie ^ 

c 



THE FOLLOWING JCL CARDS SHOULD FOLLOW JOB CARO, EACH C 

3EGINNi,N3 IN CCLUMN 1. C 

C 

// EXEC FCRTCLGV, R EGION .G0=350K C 

//FORT. SYS PR I NT DO DUMMY C 

//SYSIN DO * C 



THE SECOND OF THESE SHOULD EE REMOVED IF A PROGRAM C 

listing is DESIRED AS WELL AS OUTPUT DATA. C 

C 

C 

IMPLICIT Q£AL*8 (A-H,0-Z) 

REAL* 3 XS (25 ) , XL ( 25 ), XU (25 ) ,X( 2 ) , XDOT( 2 ) ,OELTA( 3 ) , 

*W(3) ,-^N( 3) ,B£TA(3) ,THACUT(3001 ,3) ,XDATA ( 300 1 , 5 ) , 
lA'^o(3 ),0ELAY(3) 

01 MENS I CN I0UT(3) ,IWT(3 ),ITA8(3) , IN0UT(3) 

INTEGER STE 0 ,RAMP ,PAR A 

DATA STE°,R AMO,oaRA/4HSTEO,^H°AMP ,4HPARA/ 

CCMMON THAOUT.X OAT A.T , DT ,r F , 

*NOUT, I IN ,MV ,M3, ICONT, N£C,I SKI P ,I tf,IOUT ,NQPT 
COMMON /REGl/ WtIPM 

C 

INPUT CONTROL FLAGS C 

REAO( 5,5 0 ^RUNS,NOPT,NGR APH ,NR£F , IIN,NQUT, IPRINT , I FR £0 
WRITE( 6,51 )NRUNS, IIN,NOUT 
IF{N0PT,E0.1)WRIT E( 6, 52 ) 

I=(N0PT.EQ.2) WRIT£(6,53) 

I F (NGRAPH.EO .1 )WR ITE( 6, 54) 

IF(NGRAPH.EQ.2) WR IT E( 6 ,55 ) 

IF (NGRAPH.EO, 3 ) WR ITS ( 6, 56) 

IF(NREF.£0.1 ) WR IT E( 6, 57 ) 

IF(NRE'=.5Q.2)WRITE(6,53) 

IFdPRINT.EO.l )WRITE(6,60) IFREQ 
IF( IPRI NT .E0.2 )WRIT£(6,59) 

'YO 30 I RUN=1 ,NRUNS 
5 £A0( 5,61 )T47,0T, TF 
ITF=TF/OT 

IF( ITF.GT.300i) ITF = 3C01 

IDP = ITF 

NPPLT=I0P/5 
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11 
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10 

16 

33 

32 

4 



35 

34 
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WRITE(6»62)T47,0T,TF 

r 

INPUT OPTIMIZATION DATA C 

C 

GO TO (lt2)',N0PT 

READ (5 t 50 INV,NAV,NTA,NPR, IP, IPM 
READ( 5,61 M XS( I ) ,1=1 ,NV) 

REA0(5,61 )(XU( I ),I=1,NV) 

REA0(5,61 HXL (I ),I=1,NV) 

WRI TE(6,63)NV,NAV,NTA,NPR,I P,I PN 
DO 31 1=1, NV 

WRIT = (6 ,64) I,XS( I ) 1 1 ,XU (I ) , I ,XL ( I) 

CONTINUE 

GO TO (3,4),NREF 

00 32 I=1,N0UT ^ 

input/calculate REFERENCE CURVES C 

C 

READ (5, 50)IOUT( I) ,IWT{ I) ,ITAB(I) 

W( I )=INT ( I)/100000.CO 
WRITE(6,65) I,I0U7(I),I, W (I) 

IF( I-^ABl I ) .EQ .1 ) WR I TE{ 6,76) I 
IF(ITAe(I)*E0«2)WRIT£(6,77 )I 
ITA=ITAB(I ) 

GO TO ( 5,6) , ITA 

READ (5, 61 )(XOATA( IDATA,!) , IDATA=1, ITF) 

GO TO 32 

0 EA0(5, 75 )3ET A( I ) ,DELTA( I ) , WN ( I) , AMP ( I ) ,OELAY( I ) , 
IMPUT( I ) 

-<RI TE ( 6,66) I ,3£TA ( I ) , I,DEL^A ( I ) ,I , «N (I ) ,I ,AMP( I ) , 
I, INPUT (I), I, DELAY! I ) 

T = T47 
A17=0.D0 

1 OAT A=0 
X(1 ) = 0.D0 
X (2)=0.D0 
NT = 0 

DO 33 J=1,ITF 

IF! '.GE.DSLAY! I ) .ANO.IN=»UT! I) ,£q. STEP )A17=AMP! I ) 
I F(T.GE. delay ! I ) . and. in out ! I ) . EQ«RAMP >A17=AM0! I ) 
*!T-OELAY! I ) ) 

IF!T .GE.OELAY! I ) .AND. IN PUT! I ) .EQ.PARA )A17 = AMP! I ) 
= ! !T-OELAY! I) )**2 ) 

XOOT! 1)=X!2 ) 

XDDT !2)=-2 . CO^CELTA! I ) *WN ! I 1 *X ! 2 )- ( WN ( I ) == 2) » 

X!1 )4-3ETA( I >=A17 
S = RKLD£0!2,X,XDOT,T,OT,NT) 

IF IS-1 .00 )10tll ,16 
'fiRITE! 6,6 7) I 
STOP 

XDATA! J,I ) =X!1 ) 

continue 

continue 

GO TO 12 
DO 34 1 = 1, ITF 
DO 35 J=1,N0UT 
XDATA! I , J)=0,00 
CONTINUE 
CONT I NUE 
T = T47 
I3KIP = 0 



IF SIMULATION ONLY CALL PLANT 

IF! NQPT.EC.2 ) CALL PLANT(C) 

IF! NOPT.EQ.Z) GO TO 13 

OPT I "^I 2ATIGN 

5 = 1 .DO/3 .DO 
WRITE !6,T3 ) 

CALL BOXP LX!NV, NAV, NPR,NTA, R, XS ,I?,XU ,XL,Ymn, icR ) 



C 

r 

c 

c 

r 

C 
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WRITE (6»73) 

WRITE (6, 63) 

00 36 1=1, NV 

WRIT6( 6, 69) I, XS( I ) 

36 CONTINJUE 
WRITE (6,70) YMM,IER 

PRINTEC OUTPUT 

13 GO TO (14i 15) , I PRINT 

14 WRITE(6,-^1) 

00 37 1=1, ITF ,IFREQ 

T P=DT=i«I 

WRI TE(6 ,72) TP, (XOATAd, J) , THACUT ( I , J ) , J=1 , NOUT ) 

37 CONTINUE 

PLOTTED OUTPUT 

15 GO TO (7,3,9) ,N GRAPH 

7 00 3 8 1 1 =1 , NOUT 

WR IT£(6,73) 

CALL o°LT (NPPLT , IDP, II ) 

38 CONTINUE 
GO TO 9 

8 00 39 11=1 ,NOUT 

WR I TE( 6,73) 

CALL oiC(NPPLT,IDP,II) 

39 CONTINUE 

9 WRITE(6,74> I RUN 
30 CONTINUE 

STOP 

50 FOR'^AT(8I5) 

51 FCRMAT(»0», I2,1X,»R'JN(5 ) WILL BE MftOE » , 4X, I 2, IX, 

* • INPUTS' ,4X, 12, IX, 'OUTPUTS' ) 

52 FORMATt '0* , 'OPT IM IZAT ION CALLEO FOR') 

53 FCRMAT( 'O’ ,' SIMULATION CNLY CALLEC FOR') 

54 FORMAT! ' O' ,' PRINTER PLCT(S) CALLEC FOR') 

55 FORMAL ( 'O' , 'VERSAT EC PLOTS CALLED FOR') 

56 =ORMAT( « O' , ' MO GRAPHICAL OUTatJT C.iLLED FOR') 

57 FORMAT! 'O’, 'REFERENCE CURVE CALLED FOR') 

53 FCRMAT ( * O' , ' NO REFERENCE C’JP'/E CALLED FOR') 

59 FORMAT! 'O' , 'NO PRINTED SIMULATION DATA CALLED FOR') 

60 FORMAT! 'O' , 'SIMULATION DATA PRINT FREQUENC Y = ' , IX, 12 ) 

61 FORMAT! 8 F10.5 ) 

62 FORMAT! 'O' ,' INI TI AL T I M E = ' , FI 0 . 5 , 3 X , ' STE® SIZE=',F10.5 

♦ , 3X , ' F I N'AL T I M E= ' , FI 0 . 5 ) 

63 =CRMA^( 'O' , 'NV=' , I2,3X,' NAV = ' , 12 , 3 X , ' NT A= ' , 15 , 3X , ' NPR = 
*',I5,3X, '!==', II, 3X, ' IPM=' ,11) 

64 FORMAT ( 'O' , 'XS( ' , 12, ' ) = ' , FI 0 .5 , 5 X , ' XU! ' , I 2, ’ ) = ' , FIO, 5 , 
= 5X,'XL! ’ ,12,' )=' , FI 0.5) 

65 FORMAT! 'O' , ' ID’JT! ' , 1 1 , ' ) =' , I 2 , 5X , ' W ! ' , I 1 , ' ) = ' , FI 0. 5 ) 

66 FORMAT ! ’ O' , ' BET A! ' , II , ' ) = ' , =10 . 5 , 3X , ' DEL TA ( ' , 1 1 , ' ) = ' , 

^F 10. 5,3X, • WN! ' ,11 , ' ) =' ,F10. 5,3 X, ' AMP( • , II ,' >=' ,F10 .5, 
*3X, ' INPU^( ' , II, ' ) = ' , AA ,3X, 'DELAY! ' ,11,' )=' ,F10.5) 

67 FQRMA”! 'O' ,' INTEGRATION TROUBLE-REFERENCE CURVE') 

68 FORMAT! '0 ', 'THE FREE PARAMETERS ASSOCIATED WITH THE 
= MINIMUM PERFORMANCE MEASURE ARE:') 

69 F CRMAT( • O' , 'XS! ' , 12 ,' ) =' ,F20.1 0 ) 

70 FORMAT ( 'O' , 'THE MINIMUM PERFORMANCE MEASURE IS:',2X, 
*F20.10, 5X, ' I ER=' ,2X, II ) 

71 FORMAT! ' 1' , 3X,' TI ME' , 9 X , ' XDA^ A ( 1 ) ' , 7X , ' TH AOUT ! 1 ) ' , 

* 6X, ' XDATA! 2) ' , 7 X, ' THA CUT! 2 ) ' , 6 X , ' XOA TA ( 3) ' , ?X, ' 

♦THAOUT (3 )' ) 

72 format ( 'O' , 6( Flo. 5, 5X) ,F10. 5) 

73 FQRMAT('l') 

74 FORMAT! ' O' ,' RUN NUMB £ R ' ,1X , II ) 

75 FORM AT ( 5F10 .5,A4) 

75 FORMATCO' ,' REF. CURVE! ', II. ') = ’, '"TABULATED DATA') 

77 FORMAT! 'O' ,'REF. CURVE!' ,I1,')=SECCN0 ORDER RESPONSE' ) 
END 
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SUBROUTINE PLANT (C) 

SUBROUTINE OLANT(C) SL'^ULATES THE SYSTEM 
IMPLICIT REAL-^3 (A-H,0-Z) 

REALMS G(25),P( 25 ),Z( 25),0RIVE ( 25),THA( 25) ♦0MG{ 25) ♦ 
1THA00T(25 ) t CMGOOT (25 ) , CPV I N ( 25 ) » P LAG! 25 t 25 ) » 

2 X2(25>, X2DOT(25), THACUT( 3 001 ,3 ) , 

3XDATA(3001t 3 )» C(25) , AM o ( 25 ) t DEL A Y { 2 5 ) t TAG ( 25 » 25 ) 
DIMENSION IC(25) t I0( 25 )♦ 15(25 ) ,MF (25) tIR (25)» IV(25 ) ♦ 
1IEE(25) , INPUT! 25) ,I0UT(2) 

INTEGER STEPt RAMO, PAR A 

DATA STEP* RAMP, PARA/ AhSTEP,AHRAMP,4HPARA/ 

COMMON THA0UT,XDATA,T,0T,TF, 

CNCUT, IIN,NV,M3, ICONT ,NEQ, ISKI°, ITF, IOUTtNOPT 
IF ( I SKI P-1 ) 1,5,5 
1 ISKI=> = 2 
REA0(5»2A-) N 
HI = OT 
H2 = 0.5D0=*H1 
Nil = 0 
N55 = 0 
N66 = 0 
ICK = 2*N 
WRITE (6,310) 

DO 3 1=1, N 

INPUT BLOCK DATA 

RSA0( 5i25 ) IC( I ) , ID( I ), I£{ I ) , IV( I ),G(I ), P( I ) . Z( I) 
IF! ID( I ). EO.l) Nil=Nllvl 
IF( ID( I ) .EQ.5 ) N55=N55-t-l 
IF( ID( I) . 5Q.6) N66=N66+1 

WRITE { 6,26) IC( I ) , I 0( I ) , I5( I ) , IV( I ),G( I ) ,P< I ) ,Z( I ) 

3 C CNT I iNU E 

WRI TE (6,312) 

00 316 1 = 1, II N 

INPUT FORCING -UMCTICNS 

READ! 5,311 )AMP{ I ) , DELAY! I ) , INPUT! I) , lEE ( I ) 

WRITE! 6,317) l££( I ) ,AMO(I) ,IMPUT(I) , DELAY! I ) 

6 CONTINUE 

Nil = A=N11 
N55 = 4*N55 
N66 = 4^N66 
iNEQ = N-1 

SET POINTERS 

DO 4 J=1,N 
DO 4 K=l, M 

= LAG (K, J) =0.00 

IF! IV(K).EQ. IE( J) )FLAG(K,J) =1 .00 

4 CCNTINUE 

INIT lALIZATI ON 

5 DO 6 ICLR=1 ,M 

THA( I CLP. ) = 0.D0 
THAOOT ( ICLR )=0 .00 
OMG ( ICLR) =0. DO 
OMGOQT ( ICLP )=0.00 
X2 ( ICLR)=0. DO 
X2DGT! ICLR) =0.0 0 
6 CCNTINUE 
T=0.C0 
NR2 = 1 
NR3 = 1 
NR4 = 1 
M 2 = 0 
ICONT = 0 
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IWAIT = 0 
I T = 0 
ILAST = 0 
I5LAST = 0 
I6L4ST = 0 

IF(N3PT .EO .2) GO TO 7 

ENTER C(I)=G(I) VALUES AT THIS POINT 

7 DC 9 M0RV=1,N 

DRIVE (MDRV) =0. DO 
DRVIN (MDRV )=0 .00 
00 8 ,N 

IF( IV(M).,NE.IE(M0RV) ) GO TO 315 

OR IVE (MCRV) = ORiV.E(MCRV ) +T HA { M 1*FLAG ( M,MDR V) 

315 IP ( IEE( .M) .NE. IS (MDRV) ) GO TO 8 

1= ( INPUT(M) . EQ. STEP ) GO TO 32 
IF( INPUT (M) .EQ. RAMP) GOTO 33 
IF ( INPUT(M) . EQ.PARA) GO TO 34 

32 IF( T .Gc.OELA Y(M ) )0RVIN(M0RV) =AMP(M) 

GO TO 8 

33 IF(T.Gc.OELAY(M) ) DRVIN (MDRV ) = AM P ( M ) # ( T-OEL AY (M ) ) 
GO TC 3 

34 I r (T.GE. DELAY (M) ) DRVIN (MDRV ) = AMP(M) ’i' ( ( T-0 EL AY ( M )) 
* ’^'*2 ) 

3 CONTIMUE 

DRIVE(MORV) =0RIVE( >^ORV ) +ORV I N (M r^V ) 

9 CONTIMUE 

10 M3 = M3+1 

BLOCK SIMULATION 

IP ( IWAIT.EQ.O) T = T-*-H2 
IF (I 0(M3) .EQ .1) GO TO 11 
IF ( ID( M3 ) .E0.2 ) GO TO 12 
IP ( ID(M3).EQ.3) GO TC 13 
IF ( I0( M3 ) ,eo .4) GO TO 14 
IF (ID( M3) .EQ.5 ) GO TO 15 
IF ( 10 ( M3) .EQ.6) GO TO It 
WRITE (6,23) 

STOP 

11 THA(M3) = G( M3)*0RIV£ ( M3 ) 

IWAIT = IWAlT-hl 

ILAST = ILAST+1 

IF ( ( ILAST.EQ.M 11) . ANC. ( ICONT . EG .NEC ) ) GO TO 21 
GC tq 13 

12 THAD0T(M3) = -P ( M3 ) *Th A ( M3 ) -t-G ( M3 ) =>C P I VE ( M3 ) 

3 = RKL0E2( THA, THA00T,NR2) 

IWAIT = IWAIT+1 
I F( S-1.0017,13.21 

13 0MG00T(M3) = -? ( M 3 ) =^QMG( M3 ) +0 ( M3 ) *0 RI VE ( M3 ) 

5 = RKLDE3 ( CMG, OMGDOT ,NR3) 

THA(M3) = ( Z(M3)-P( M3) )»GMG(M3 ),.G(M3)*0RIVS(M3 ) 

IWAIT = IWAIT + 1 
IF(S-1.00)17,ia,21 

14 V = :: PLX( P,Z,G,0RIVE,X2,0MG,NR4) 

THA(M3) = 0MG(M3) 

IWAIT = lWAIT+1 
IF(V-l.DC) 17, 13,21 

15 THA(M3) = G(M3)»ORIV£(m2 ) 

1= (THA(M3) .GT.P( M3 ) ) TH A( M3 ) =P ( M3) 

IF (THA(M3) .LT.Z(M3) ) THA ( M3 ) = Z ( M 3 ) 

IWAIT = IWAIT-t-1 
I BLAST = I 5LAST+1 

IF (( I5LAST.EQ.N55) .ANO.(ICONT.EQ.NEQ) ) GO tq 21 
GO TO 18 

16 THA(M3) = G(M3) ^ORI VE( M3) 

I F( ( DABS (T HA( M3 ) ) ).LT.P(M3) ) THA ( M 3 ) =0 .0 0 
IWAIT = IWAIT4-1 
I BLAST = I BLAST vl 

IF ( ( I6LAST.EQ.N66) .ANC.( ICONT.EQ .NEC) ) GO TO 21 
GO TO 1 a 
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13 

19 



20 

21 

38 
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23 
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25 
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27 
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29 

30 



310 

311 

312 
217 



11 

2 

3 
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WRIT? (6,29) 

STOP 

ICON! = ICONT+l 

IF (ICONT-N) 10,19,20 

■NR2 = NR2+1 

NR3 = NR3-H 

NR4 = NR 4+1 

M3 = 0 

ICQNT = 0 

IF { IWAIT.EO.ICK) IWAIT=0 
GO TO 7 
WRITE (6,30) 

STOP 

STORE OUTPUT DATA 

IT = IT+1 
DC 33 I=1,NCUT 

THAOUTUT,!) =THA(I0UT(I) ) ^ 

CONTINUE 

IF (T-TF) 22,22,23 

NR2 = 1 

NR3 = 1 

NR4 = 1 

M3 = 0 

ICON"^ = 0 

IWAIT = 0 

ILAST = 0 

I5LAST = 0 

I 6LAST = 0 

GO TO 7 

RETURN 

FORMAT! 12) 

FORMA* (412, 2X,3F1 0.5 ) 

= CRMAT (/,5H 8LK ,12,I1,2H = ,2 12 , 6X, 3 =20 .7 ) 

FORMAT (//, 2X, • THETA OUT IS A ( * ,I 2 , ' J • ) 

FORMAT (40H EON SWITCH CONTROL DID NOT WORK ***) 

FORMAT (20H INTEGRATION TROUBLE) 

FGRmaT(58X, 'ERROR IN INTERGA-^I ON. » , 2X, 

I'ATTEMPTED TO INTEG. '^CR£ THAN N-EGTN* ) 
format (• O' , 40X, » THE SYSTEM BLOCKS ARE:') 

FORMAT (2F10.5,A4, 12 ) 

FORMAT (• C ,40X, • THE SYSTEM FORCING FUNCTIONS ARE:') 
FORM A’ ( *0 ' , 'CRV IN ( ' , 12, ' )=’ ,F10. 5,A4,F10. 5) 

E NO 



FUNCTION RKLOEQ (N ,X , X DOT , T , DT , NT ) 

CALCULATES RESPONSE OF SECOND 
ORDER REFERENCE CURVE 

implicit REAL*8 (A-H,0-Z) 

REALMS X (2) ,X DOT ( 2 ) ,0 ( 25 ) 

NT = NT + 1 

GO TO ( 1,2, 3,4) ,NT 

HI = DT 

H2 = HlTC.500 

H3 = HI *2.000 

H6 = Hl/6.000 

00 11 J = 1,N 

Q(J) = 0.00 

A = 0.500 

T = T + H2 

GO TO 5 

A = 0,2928932133134525 
GO TO 5 

A = 1 .7071 06731 1365475 
T = T + H2 
GO TO 5 
DC 41 I = 1,N 

X(I) = X(I) + H6^X00T(I) - Q(I)/3.00 
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NT = 0 
RKLDEO = 2. 

GC TO 6 

5D0 51L=1,N 

X(L) = X(L) + DT*XDOT( L)-0(L ) ) 

51 Q(U = H3*A*XDOT(L) + (l.CO - 3.00=*A)^Q(U 
RKLD6Q = 1. 

6 RETURN 
END 



FUNCTION RKLDE2 (X,XD0T,NR2J 
IMPLICIT R£AL*3 (A-HtO-Z> 

REALMS THAQUT(3001,3) » XOATA ( 3001,3) 

REAL*3 X(^), X00T(4), G(25) 

DIMENSION I0UT(3) 

common THA0UT,XDATA, T,DT, TF, 

CNCUT,II N,NV,M3,ICaNT,NSC,ISf<IP, ITF, lOUT 
GO TO ( 1,2, 3,4) , NR2 

1 HI = OT 

H2 = Hi»0.5D0 
H3 = Hl*2.000 
H6 = H1/6.000 
Q(M3) = 0.00 
A = 0. 5D0 
GG TO 5 

2 A = 0.2528932138134525 
GO TO 5 

3 A = 1 .7071067811965475 
GO TO 5 

4 X(M3) = X( M3)+H5^XOOT(M2)-Q(M3) /3.0C 
RKL0E2 = 1. 

IF { ICONT.E G.NEQ) RKLDE2=2. 

GO TO 6 

5 X(M3) = X(M3) -t-A=i‘ ( CT*XOCT (M3)-G( M3 ) ) 

Q(M3) = H5=«A*XD0T( M3)+(1.00-3.00*A>*Q( M3) 
PKLDE2 = 1. 

6 RETURN 
END 



FLNCTION RKLDE3 (X,X0CT,NR3) 

IMPLIC IT peaL*3 (A-H,0-2) 

REAL*8 THAOUT (3301 ,3 ) , XCAT 4 ( 300 1 , 3 ) 

REAL’^S X(4), XOOT(4), Q(25) 

0IM6NS ICN I OJT ( 3) 

CCMMCN THAOUT, XOATA, T ,CT ,T=, 

CNOUT, I IN,NV ,M3, ICONT,NEQ,I SKI P , ITF, lOUT 
GC TO (1,2, 3,4W NR3 

1 HI = OT 

H2 = Hl*0.500 
H3 = HI *2 .000 
H 6 = HI /6. 000 
Q (M3 ) = 0.00 
A = 0.500 
GO TO 5 

2 A = 0.2923932183134525 
GC TO 5 

3 A = 1.7071067311365475 
GC TO 5 

4 X(M3) = X( M3) +H6*XD0T (M3 )-G(M3 ) /3 .CO 

RKL0E3 = 1. 

IF ( IC0^7,EC.NEQ) RKL0£3=2. 

GC TO 6 

5X('^3) = X(M3) +A«( 0T*XD0T(M3)-Q( M3) ) 

Q(M3) = H3>^A»X00T (M3) +(1 .00-3 .OO^A )»0(M3) 
RKLDE3 = 1. 

6 RETURN 
END 



FUNCTION CCPLX ( P , Z , C- , CR IV E , X2 , 0 MG, NR4 ) 
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implicit realms (A-HtQ-Z) 

REALMS THAOUTOOOl »3 ) » XC AT a ( 3 001 , 3 ) rCMGDGT(l) 

REAL-1'8 0MG{ 1 ) ♦ P ( 1 > , Z { 1 ) t G ( 1 ) ,0 R I VE ( I) ,X 2 ( 2 5 ) ♦ X2D CT ( 25 ) 

DIMENSION I0UT(3J 

COMMON THAOUT,XDATA»TtOT,TF, 

CNOUT, IINtMV,M3» ICGNT,N£q» I SK I P , I TF , lOUT 
0MGD0T(iM3) = X2 (M3) 

SS = RKL0E3( OMG ,0MG0QT,NR4) 

X2D0T{.M3) = -2.=*P(M3)*Z(M3)*X2(M2)-ZIM3)=»^2*0MG(M3) + 
LG(M3)*0RIVE( M3 ) 

SSS = RKLDE4(X2»X200T,NR4) 

CCPLX = SSS 

RETURN 

END 



FUNCTION RKLDE4 (X»X0CT»NR4) 

IMPLICIT REALMS (A-H,0-Z) 

R£AL*8 THAOUT (3001,3 ) , X OAT A { 300 1 , 2 ) 

REALMS X(l) , XDOTd), CC(251 
DIMENSION IQUT{3) 

CC'^MGN THAOUT ,X DATA, T ,DT,TF, 

CNOUT, IIN,NV ,M3,IC0NT,N£g,ISKIP,ITF,I0UT 
GO TO ( 1,2, 3,4) , NR4 

1 HI = DT 

H2 = H1#0.5D0 
H3 = HI ^2 .000 
H6 = Hl/6.000 
0C(M3) = O.DO 
A = 0 .5D0 
GC TO 5 

2 A = 0. 2S28932138134525 
GC TO 5 

3 A = 1,7071067811365475 
GO TO 5 

4 X(M3) = X( M3 )-Hi6tXOOT (M3 )-QC(M3 )/ 3.00 
RKLDE4 = 1. 

IF { ICONT.EQ.NEQ) RKLD£4=2. 

GC TO 6 

5 X(M3) = X(M3)+A=s{DT*XD0T(M3)-qC( M3) ) 

QC(M3) = H3»A^XD0T(M3 )+( 1.00-3. D0*A)*QC(M3) 

RKLD64 = 1 . 

6 RETURN 
END 



FUNCTION XE(C) 

EVALUATES IMO LICIT CCNSTRAINTS FCR SCXPLX 

REAL*S C(25) 

KE=0 
R ETUR.N 
ENID 



FUNCTION FE(C) 



THE PUNCTION MINIMIZED BY BCXPLX 



i mplicit REAL»8 ( A-H.O-Z ) 

REALMS THA0UT(3Q01,3) ,XDATa(3 00: ,3) 
COMMON THAOUT, XDATA, T ,OT,^F, 

*.NOUT,I I N,NV ,M3, ICQNT , NSC. ISKIP, IT-. 
COMMON /PEGl/ W,IPM 
DC 10 I=1,NCUT 
THAOUTd ,I ) =0.00 
10 CONTINUE 
OIFF=O.CO 
PI =0.0 0 
CALL °LANT(C) 

T =0T 



,C (25 ) ,'w (3 ) 
I'^UT 
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DC 11 1=1, ITF 
00 12 J=l tNOUT 

0IFF=XD4TA( I, J )-THA0UT( I , J ) 

GO TQ(l ,2,3,^ ), IPM 

1 PI =PI+-W (J)* ( 0IFF*2) 

GO TO 12 

2 PI = PH-W (J )* ( 0IFF**2 
GO TO 12 

3 PI = PI+W(J)*OABS(OIFF) 

GO TO 12 

A PI =PI+W( J)*OA3S(OIFF )*T 

12 CONTINUE 
T=T+OT 
11 CONTINUE 
FE=P I 
RETURN 
END 

✓ 

SUBROUTINE PPLT ( N PPLT , JCP, 1 1 ) 

AUTOSCALES AND CALLS FOR PRINTR ©LOTS 

IMPLI Cl T REALMS (A-H,C-2) 

REAL^A XX( 900 ), YY (900 I , VW( 900) 

R£AL*A TX(A),TY(A) 

REAL *3 THA0UT(300 1,3) , XDATA ( 3 001 , 3 ) 

Olf^ENSlCN I GUT (3) 

common THA0UT,X0ATA,T ,DT ,TF, 

CNOUT, 1 1 N ,NV ,M3 ,ICONT, NEC, I SKI ?, ITF, I OUT 
OC 1 1=1,900 
XX( I ) = 0, 

YY( I ) = 0. 

1 WW(I) = Q. 

T=0. 00 

TSTEP=5.00*0T 

J = 0 

3 IGX=O.DO 

3 IGY = 0 .DO 

SMLX=0.00 

SMLY=O.DO 

DC 2 I=1,IDP,5 

J = J + 1 

XX( J ) = T 

YY( J)=XOATA(I,II ) 

WW( J) =THAQUT(I ,11) 

X = T 

XO = YY (J) 

TH = WW( J) 

YMAX=0MAX1 (XDtTH) 

YMIN=OMIM (XD,TH) 

XyAX=DMAXl (BIGX ,X) 

IF (B IGX.LT.X) 3IGX=X 
IF (BIGY.LT.YMAX) 3IGY = Y'<AX 
IF ( SML Y.GT.YMIN) SMLY=YMIN 

T = t+tsTEP 

2 CONTINUE 
TX(1) = 0. 

TX(2 ) = 0. 

TX(3) = SMLX 

TX(A) = BIGX 
TY(1 ) = 3IGY 
TY(2) = SMLY 

TY(3I = 0. 

TY(A) = 0. 

WRITE (6,3) BIGX, 3IGY ,SMLY 
CALL °LOTo (TX,TY,A, 1 ) 

CALL PLCTO (XX,WW,NPPLT ,2) 

CALL PLOT? (XX,YY,NPPLT,3) 

W RITE (6 ,A ) 

I F( IDP. EQ.A500) WRITE(6,5) 

PETURN 
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3 FORMAT ( 2X, ' BIGX= • t El 5 . 7 ,2 X, • B I GY = '♦E15.7t2Xt 
l'SMLY= NE15.7) 

4 PQRMAK //, 2 X,» SYSTEM RESPONSE FOR PRCBLEV ’) 

5 format ( //,2X».' STOP AT 900 GRAPH POINTS’) 

EKO 



SUBROUTINE PIC (NOPLT t IDP, 1 1 ) 

AUTOSCALES AND CALLS FOR VERSATEC PLOTS 

IMPLICIT REAL=^3 (A-H,C-Z) 

PEAL *4 XX(900),YY(900),V»W(900) 

R£AL*4TX(4) tTYI 4) 

REALMS THA0UT(3001,3) , XOATA ( 3 001 » 3 ) t TI TLE ( 12 ) 

REAL *3LABC/ • • / 

REAL *4LABA/' A ’/,LA8C/‘ 0 •/ 

DIMENSION I0UT(3) 

CCMMON THAOUTtX OATAtT ,DT,TF» ' 

*NOUTtIIN,NV ,M3,ICaNT, NE C 1 1 SKI = t IT F , I OUT 
READ (5,2) TITLE 
T=0 .00 

TSTEP =5.D0*0T 
J = 0 

3IGX = 0.00 
3IGY = o.no 
SMLX = 0.00 
Sf'LY = 0.00 
00 1 I =1 ,IOP ,5 
J = J+1 
XX(J) = T 
YY( J ) = X0ATA (I tl I ) 

WW(J )=^HAQUT( I, II ) 

X = T 

XD = YY( J) 

Th = WW(J) 

YMAX = CMAXKXOtTh) 

YMIN = OMINl ( XD ,TH) 

X.MAX = CMAXK BIGXtX > 

IF (BIGX.LT.X) 3IGX=X 
IF (BIGY.LT.YMAX) 3IGY=YMAX 
IF (SMLY .GT.YMIN ) SMLY=YMIN 
T = T+TSTE? 

1 continue 

^X(l) = 0. 

TX(2) = 0. 

*X(3) = SMLX 
TX(4) = aiGX 

TY(1) = BIGY 
TY(2) = SMLY 
TY(3 ) = 0. 

TY(4) = C. 

CALL DRAW ( 4 , TX , T Y , 1 , 1 , L ABC , T I TLE , C , C, C, C , C, 0, 3 , 3 , 0 , L 
CALL DRAW (NPPLT,XX,WW ,2 ,0,LAaA,TrLE,0,0,0,0,OtO,3,8 
^0 t L ) 

CALL DRAW ( NPPLT, XX, Y Y , 3 , 0 , L A 3 0 , T ITL E , C , C, 0 , 0, 0 , 0, 3 , 3 
’*0,L) 

IF (L.NE.O) WRITE (6,3) L 
RETURN 

2 =0RMAT (6A3) 

3 FOR'^A^ (//,' GRAPH NOT COMOLETEO. OUTPUT CODE = ',12 
END 



SUBROUTINE B0X3LX ( NV , NA V , NPR , NT Z , R Z, XS , I P , 3U , 3 L , YMN , 
♦I ER) 

It * series; :*Sc * :^c3c * :S -ft 

SUBROUTINE 80XPLX (CATEGORY hO ) 

PUROQSE 
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30XPLX IS A SUBROUTINE USED TO SOLVE THE OROBLEH 
OF LOCATING A MINIMUM (OR MAXIMUM). OF AN ARBITRARY 
OBJECTIVE FUNCTION SUBJECT TO ARBITRARY E>PLICIT 
AND/OR IMPLICIT CONSTRAINTS BY THE COMPLEX METHOD 
OF M.j. BCX. EXPLICIT CONSTRAIN'^S ARE DEFINED AS 
UPPER AND LOWER BOUNDS ON THE INDEPENDENT VARIABL- 
ES. IMPLICIT CONSTRAINTS MAY BE ARBITRARY P'JNCTIONS 
OF THE VARIABLES. TWO FUNCTION SUBPROGRAMS TC 
EVALUATE THE OBJECTIVE FUNCTION AMO IMPLICIT CGN- 
STRAINTSt RESPECTI VELY f MUST BE SUPPLIED BY THE 
USER (SEE EXAMPLE BELOW). 30XPLX ALSO HAS THE OPT- 
ION TO PERFORM INTEGER PROGRAMMING, WHERE THE VAL- 
UES OF THE INDEPENDENT VARIABLES ARE RESTRICTED TC 
INTEGERS. 

USAGE 

CALL BGXPLX (NV ,N AV , NPR ,NT- A . R . X S t IP, XU, XL t YMN , lER ) 
DESCRIPTION OF PARAMETERS 

NV AN integer INPUT DEFINING THE NUMBER OF INDE- 
PENDENT VARIABLES OF THE OBJECTIVE FUNCTION 
TO BE MINIMIZED. NOTE: MAXIMUM NV + NAV IS 
PPESENTLY 50. MAXIMUM NV IS 25. IF THESE 
LIMITS MUST BE EXCEEDED, PUNCH A SOURCE DECK 
IN THE USUAL MANNER, AND CHANGE THE DIMEN- 
SION STATEMENTS. 

NAV AN INTEGER INPUT DEFINING the NUMBER OF AUXI- 
LIARY VARIABLES THE USER WISHES TO DEFINE FOR 
HIS OWN CONVENIENCE. TYPICALLY HE MAY WISH 
TO DEFINE THE VALUE OF EACH IMPLICIT CON- 
STRAINT FUNCTION AS AN AUXILIARY VARIABLE. IF 
THIS IS DONE, THE OPTIONAL OUTPUT FEATURE 
OF BOXPLX CAN BE USED TO OBSERVE THE VALUES 
OF THOSE CONSTRAINTS AS THE SOLUTION PRO- 
GRESSES. AUXILIARY VARIABLES, IF USED, SHOULD 
BE EVALUATED IN FUNCTION KE (DEFINED BELOW). 
NAV MAY 3c ZERO. 

NPR INPUT INTEGER CONTROLLING THE FREQUENCY OF 
OUTPUT DESIRED FCR DIAGNOSTIC PURPOSES. 1 - 
NPR .LE. 0, NO OUTPUT WILL BE PRODUCED BY 
BOXPLX. OTHERWISE, THE CURRENT COMPLEX OF 
K = 2’!‘NV VERTICES AND THEIR CEN'ROID WILL BE 
OUTPUT AFTER EACH N°R PERMISSIBLE “TRIALS. THE 
NUMBER OF TOTAL TRIALS, NUMBER OF FEASIBLE 
TRIALS, NUMBER OF FUNCTION EVALUATICNS AND 
NUMBER DF IMPLICIT CONSTRAINT EVALJA'^iONS ARE 
I NCLUOED IN THE OUTPUT . 

ADDITIONALLY, (WHEN N®R .GT. 0) THE SAME IN- 
FORMATION WILL BE OUTPUT; 

1) IF THE INITIAL POINT IS NOT FEASIBLE. 

2) AFTER THE FIRST COMPLETE COMPLEX IS 
GENERATED. 

3) IF A FEASIBLE VERTEX CANNOT BE FOUND AT 
SOME TRIAL, 

4) IF THE OBJECTIVE VALUE CF A VERTEX CAN- 
NOT BE MAGE NQ-LONGER-WORST . 

5) IF THE LIMIT ON TRIALS (NTA) IS REACHED 
AMD , 

6) when the OBJECTIVE FUNCTION HAS BEEN UN- 
CHANGED FOR 2«NV TRIALS, INDICATING A 
LOCAL MINIMUM HAS BEEN FOUND. 

1= THE USER WISHES TO TRACE ^he PRCGRES3 OF 
A SOLUTION, A CHOICE OF NPR= 25, 5C OR 100 IS 
RECCMMENDEO, 



NTA INTEGER INPUT OF LIMIT ON THE NUMBER OF 
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TRIALS ALLOWED IN THE CALCULATION. IF THE 
USER INPUTS NTA .L = . 0, A CcFAULT VALUE OF 
2000 IS USED. WHEN THIS LIMIT IS REACHED 
CONTROL RETURNS TO THE CALLING PROGRAM WITH 
THE EEST ATTAINED OBJECTIVE FUNCTION VALUE IN 
y^N, AND THE BEST ATTAINED SOLUTION OCINT IN 

xs. 

R A REAL NUMBER INPUT TO DEFINE THE FIRST RAN- 
DOM NUMBER USED IN DEVELOPING THE INI- 
TIAL COMPLEX OF 2*NV VERTICES. 

(0. GT. R .LT. 1.) IF R IS NOT WITHIN THESE 
BOUNCSt IT WILL BE REPLACED BY 1./3. . 

XS INPUT REAL ARRAY DIMENSIONED AT LEAST NV+NAV. 
THE FIRST NV MUST CONTAIN A FEASIBLE ORIGIN 
FOR STARTING THE CALCULATION. THE LAST NAV 
NEED NOT BE INITIALIZED. UPON RETURN FROM 
BOXPLX, THE FIRST NV ELEMENTS OF THE ARRAY 
CONTAIN THE COORDINATES OF THE MINIMUM OB- 
JECTIVE FUNCTION, AMO THE REMAINING NAV 
(NAV ,GE. 0) CONTAIN THE VALUES OF Tj-e 
CORRESPONDING AUXILIARY VARIABLES. 

IP INTEGER INPUT FOR OPTIONAL INTEGER PRO- 
GRAMMING. IF IP=1, THE VALUES OF THE INDE- 
PENDENT VARIABLES WILL BE REPLACED WITH 
INTEGER VALUES (STILL STCREO AS REAL*4I. 

XU A REAL ARRAY DIMENSIONED AT LEAST MV INPUT- 
TING THE UPPER BOUND ON EACH INDEPENDENT 
VARIABLE, (EACH EXPLICIT CONSTRAINT). INPUT 
VALUES ARE SLIGHTLY ALTERED BY BOXPLX. 

XL A REAL ARRAY OIMENSIOMcD AT LEAST NV INPUT- 
TING THE LOWER BCUNO ON EACH IN0EP5N0ENT 
VARIABLE, (EACH EXPLICIT CCNSTRAINT). NOTE: 
FOR BOTH XU AND XL CHOOSE REASONABLE VALUES 
IF NONE ARE GIVEN, NOT VALLES WHICH ARE MAG- 
NITUDES ABOVE OR BELOW THE EX0£CT5D SOLUTION. 
INPUT VALUES ARE SLIGHTLY ALTERED BY BOXPLX. 

YMN THIS OUTPUT IS THE VALUE (REAL^A) OF th£ OB- 
JECTIVE FUNCTION, CORRESPONDING TO THE SOLU- 
TION POINT OUTPUT IN XS. 

I6R INTEGER ERROR RETURN. TO EE INTERROGATrC 

UPON RETURN FROM BOXPLX. lER WILL 35 ONE OF 
THE FOLLOWING: 

=-l CANNOT FIND FEASIBLE VERTEX CR FEAS- 
IBLE centroid AT THE START OR A RE- 
START (SEE 'MET HOC BELOW). 

=0 FUNCTION VALUE UNCHANGED FOR *N' 

TRIALS. (WHERE N = 6’*NV+10) THIS IS 
The NORMAL RETURN PARAMETER. 

= 1 CANNOT DE'yELOP FEASIBLE VERTEX. 

=2 CANNOT DEVELOP A NO-LCNGER-WQRST VER- 
TEX. 

= 3 LIMIT ON TRIALS REACHED. (NTA EX- 
CEEDED) 

NO^E: VALID RESULTS MAY BE RETURNED IN 

ANY OF THE ABOVE CASES. 

example of USAGE 

THIS EXAMPLE MINIMIZES THE OBJECTIVE FUNCTION 

SHOWN IN THE EXTERNAL FUNCTION F£(X). THERE ARE 

TWO INDEPENDENT VARIABLES X (1 ) S X(2), AND TWO IM- 

PLi: I T CONSTRAI NT FUNCTIONS X(3) i X(4) WHICH ARE 

EVALUATED AS AUXILIARY VARIABLES (SEE EXTERNAL 

PIJNCTI ON KE (X ) ) . 
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DIMENSION XS (4) ,XU( 2) ,XL( 21 



STARTING GUESS 
XS(1) = 1.0 
XS(2) = 0.5 
UPPER LIMITS 
XU(1) = 6.0 
XU(2 ) = 6. 0 
LOWER LIMITS 
XL(l > = 0.0 
XL(2) = 0.0 



R = 9./15. 
NTA = 5000 
NPR = 50 
NAV = 2 
. NV = 2 
lo = 0 



CALL BOXPLX ( N V , NA V , N PR tNTA t R t XS , I P » XU , XL t YMN , I £R ) 
WRITE(6 ,1 ) ((XS( 1 ) t 1= If A) ,YMNf I ER ) 

1 F0RMAT( //// f' THE POINT IS LOCATED AT (XS(H=) 
=^'A(E13.7,5X» , //, • AND THE FUNCTION VALUE IS », 
*613.7 , • lER = ',15) 

STOP 

END 



PUNCTICN K6 IX) 

EVALUATE CONSTRAINTS. SET K£=0 IF NC IMPLICIT 
CONSTRAINT IS VIOLATED, CR SET KE=1 IF ANY IMPLICIT 
CONSTRAINT IS VIOLATED. 

DIMENSION X(4) 

XI = X ( 1 ) 

X2 = X(2) 

KE = 0 

X(3 > = XI + 1.732051*X2 

IF (X(3) .LT. 0. .OR. X(3) . GT . 6.) GO TO 1 
X(L) = Xl/1. 732051 -X2 
IF (X (4) .GE. 0. ) RETURN 

1 KE = 1 - 

RETURN 
END 



FUNCTION FE(X) 

DIMENSION X(4) 

THIS IS THE OBJECTIVE FUNCTION. 

FE= “(X(2)’»*3 * (9.-( X( l)-3. )*=*2 )/(46.76533) ) 

RETUR N 
END 

METHOD 

THE COMPLEX METHOD IS AN EXTENSION AND ADAPTION OF 
’HE simplex method CF LINEAR ^RPORAMMING. START- 
ING WITH ANY ONE FEASIBLE POINT in N-DIMENSION 
SPACE A ’’COMPLEX" OF 2*N VERTICES IS CCNSTRUC’ED 
BY SELEC'^’ING RANDOM =0 IMTS WITHIN th£ FEASIBLE 
REGION. FOR THIS PURPOSE N CGCFCINATES AR£ FIRST 
RANDOMLY CHOSEN WITHIN THE S=>ACE BOUNDED BY EX- 
PLICIT CONSTRAINTS. ThIS DEFINES A TRIAL INITIAL 
VERTEX. IT IS THEN CHECKED FOR POSSIBLE VIOLATION 
CF IMPLICIT CONSTRAINTS. 1= ONE OR more aRE VIO- 
LATED, THE TRIAL INITIAL VERTEX IS DISPLACED HALF 
OF ITS DISTANCE FROM THE CSNTRQIO Qp PREVIOUSLY 
SELECTED INITIAL VERTICES. IF NECESSARY THIS CIS- 



69 



OOOOOOC ^oc^oooooooooooooooooc^ooooonoooo^'^^^ooooooooooooooooooooooooooooooooo 






p 



I 



i 



♦ 




000o0000r>000r^0000o0000000o0000c^000o^^000o0o0000c^00000000000o000«^0000(^0000 



PLACEMENT PROCESS IS REPEATED UNTIL THE VERTEX HAS 
BECOME FEASIBLE. IF THIS FAILS TC HAPPEN AF^ER 
5#N+10 OISPLACEMENTSt THE SOLUTION IS ABANDONED. 
AP'^ER EACH VERTEX IS ADDED TO THE COMPLEX, THE 
CURRENT CENTROID IS CHECKED FOR FEASIBILITY. IF 
IT IS INFEASIBLE, THE LAST TRIAL VERTEX IS ABAND- 
ONED AND AN EFFORT TC GENERATE AN ALTERNATIVE 
TRIAL VERTEX IS MADE. IF 5=J'N+1C VERTICES ARE 
ABANDONED CONSECUTIVELY, THE SOLUTION IS TER- 
MINATED. 

IF AN INITIAL COMPLEX IS ESTABLISHED, THE BASIC 
COMPUTATION LOOP IS INITIATED. THESE INSTRUCTIONS 
FIND THE CURRENT WORST VERTEX, THAT IS, THE VERTEX 
WITH THE LARGEST CORRESPONDING VALUE FOR THE OB- 
JECTIVE FUNCTION, AND REPLACE THAT VERTEX BY ITS 
OVER-REFLECTION THROUGH THE CENTROID OF ALL OTHER 
VERTICES. (IF THE VERTEX TQ BE REPLACED IS CON- 
SIDERED AS A VECTOR IN M-SPACE, ITS OVER- 
REFLECTION IS OPPOSITE IN DIRECTION, INCREASED IN 
LENGTH BY THE FACTOR 1.3, AND CCLLINEAR WITH THE 
REPLACED VERTEX AND CENTROID OF ALL OTHER 
VERTICES. ) 

when an over-reflection IS NOT FEASIBLE OR REMAINS 
WORST, IT 1$ CONSIDERED NQT-P ER V I SS I 2L E AND IS 
DISPLACED halfway TOWARD thE CENTROID. AFTER FOUR 
SUCH ATTEMPTS ARE MACE UNS JCCESSFULLY , EVERY FIFTH 
ATTEMPT IS MADE BY REFLECTING THE OFFENDING VERTEX 
THROUGH THE PRESENT BEST VERTEX, INSTEAD OF 
THROUGH the centroid. IF 5*N+10 DISPLACEMENTS AN C 
OVER-REFLECTI ONS OCCUR WITHOUT A SUCCESSFUL 
(»ER»^ISSI3L£) RESULT, THE CURREN’T BEST VERTEX IS 
TAKEN AS AN INITIAL FEASIBLE POINT FOR A RESTART 
RUN OF THE COMPLETE PROCESS. RESTARTING IS ALSO 
undertaken when 6*NV+10 CONSECUTIVE TRIALS HAVE 
BEEN MADE WITH NO SIGM=ICAM“ CHANGE IN THE VALUE 
OF THE OBJECTIVE FUNCTION. IN ALL CASES, RESTART- 
ING IS INHIBITED IF THE LAST RESTART CIO NOT PRO- 
DUCE A SIGNIFICANT IMPROVEMENT IN THE MINIMUM 
ATTAINED. 

IT IS RECOMMENDED THAT THE USER READ THE REFER- 
ENCE FOR FURTHER USEFUL INFORMATION. IT SHOULD BE 
NOTED THAT THE ALGORITHM DEFINED THERE HAS BEEN 
ALTERED TO FIND THE CONSTRAINED 'MINIMUM, RATHER 
THAN THE MAXIMUM. 



REMARKS 

THE INTEGER PROGRAMMING OPTION WAS ADDED TO THIS 
PROGRAM AS SUGGESTED IN REFEPcNCE (2). A MIXED 
riTEGER/CONTlNUOUS VARIABLE VERSION CF BCXPLX 
WOULD BE EASY TO CREATE BY DECLARING "IP" TO BE AN 
ARRAY OF NV CONTROL VARIABLES WHERE IP(I)=1 WOULD 
INDICATE that THE I-TH VARIABLE IS TO St CONFINED 
TO integer values. EACH STATEMENT OF THE FOR*^ 

'1= (IP .EQ. 1) • ETC. WOULD ^HEN NEED tq bE AL- 
TERED TO 'IF (IP(I> .EO. 1)' ETC., WHERE ^HE SUB- 
SCRIPT IS APORO°R lAT ELY CHOSEN. NORMALLY, XU AND 
XL VALUES ARE ALTERED TO BE AM E = SILON 'WITHIN' 
ACTUAL values DECLARED BY THE USER. THIS ADJUST- 
MENT IS NOT MADE WHEN I°=l. 

MOTE: NO NQN-LINEAR 'PROGRAMMING ALGORITHM CAN 

GUARANTEE THAT THE ANSWER =OUMD IS THE GLOBAL 
MINIMUM, RATHER THAN JUST A LOCAL 'MINIMUM. HOW- 
EVER, ACCORDING TO REF. 2, THE COMPLEX ’^ETHOD HAS 
AN advantage IN THAT IT TENDS tp cj.sjg GLOBAL 

MINIMUM MOPE FREQUENTLY THAN MANY OTHER NCN- 
LINEAR programming AL COR ITH'^S . 
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IT SHCULD BE iNOTEO Th>lT THE AUXILIARY VARIABLE 
FEATURE CAN ALSO 3E UiEO TD DEAL WITH RRGBLEMS 
CONTAINING EQUALITY CONSTRAINTS. ANY EOUALTIY 
CONSTRAINT IMPLIES THAT A GIVEN VARIABLE IS NOT 
TRULY INDEPENDENT. THEREFORE t IN GENERAL » ONE 
VARIABLE INVOLVED IN AN EQUALITY CONSTRAINT CAN 
BE RENUMBERED FROM THE SET OF NV INDEPENDENT 
VARIABLES AND ADDED TO THE SET OF NAV AUXILIARY 
VARIABLES. THIS USUALLY INVOLVES RENUMBERING 
THE INDEPENDENT VARIABLES OF ^hE GIVEN PROBLEM. 

SLBROUTINES AND FUNCTIONS REQUIRED 

SUBROUTINE ‘BOUT' ANC FUNCTION 'Fev* ARE INTEGRAL 
PARTS OF THE SOXPLX PACKAGE. 

TWO FUNCTIONS MUS’^ BE SUPPLIED BY THE USER. THE 
FIRST, KE(X), IS USED TO EVALUATE ^HE IMPLICIT 
CONSTRAINTS. SET K E=0 AT THE BEGINNING OF THE 
FUNCTION, THEN EVALUATE THE IMOLICIT CONSTRAINTS. 
IN THE example ABOVE, THE FIRST CONSTRAINT, X(2), 
MUST BE WITHIN THE RANGE (0. .L E . X(3) .LE. 6). 
THE SECONC CONSTRAINT X(^), MUST BE .GS. 0. . IF 

either CONSTRAINT IS NOT WITHIN THESE 3CUN0S, 
CONTROL IS TRANSFERRED TO STATEMENT 1, ANC K5 IS 
SET TO n ” AND CONTROL IS RETURNED TO BCXPLX. 

THE SECOND FUNCTION TH6 USER '^UST OROV IDE EVALU- 
ATES THE OBJECTIVE FUNCTION. IT IS CALLED F£(X) 
AS SHOWN IN THE EXAMPLE ABOVE, AND FE MUST 3£ SET 
TO THE VALUE OF THE CBJECTIVE PUNCYION CORRES- 
PONDING TO CURRENT VALUES OF THE NV INDEPENDENT 
VARIABLES IN ARRAY ’X*. 
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CORRECTED 3/1976 



IMPLICIT REAL^ (A-H,0-Z) 

REAL* 3 V(50 ,50) ,FUN( 30 ) ,SUM(25 ) ,CEN (25) ,XS (NV ) , BU( NV ) , 
*BL( NV) 

KV = 5 
£P = l.D-6 
N^A = 2000 

IF (NTZ.GT.O) NTA = NTZ 
R = R Z 

IF(R.L£.0.C0 .or .R .GE .1 .CO )R=1 .DO/ 2 .DO 
N VT = NV+NAV 

TOTAL VARS, EXPLICIT ^LUS IMOliCIT 

NT = 0 

CURP ENT TRIAL NO. 

NPT = 0 

C CURRENT NO. OF PERMISSIBLE TRIALS 

NTFS = 0 
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CURRENT MG, OF TIMES F HAS BEEN ALMOST UNCHANGED 

CHECK FEASIBILITY OF START POINT 

DC A 1=1, NV 
VT = XS(I ) 

IF ( SL ( I) .LE.VT ) GO TO 1 
II = -I 
VT = 3L( n 
GO TO 2 

1 IF (BUI I) .GE.VT ) GO TC 3 
II = I 

VT = BU( n 

2 IF (NPR.GT.O) WRITE (6,49) II 

3 V(I ,1) = VT 

C EN ( I ) = VT 

IF ( IP. EQ.l ) GO TO 4 

BL(I )=BL(I) +OMAXK EP , £P*DABS( BL ( I) ) ) 

BUI I ) =BU( I )-0-MAXl ( £P, EP’^CABS I BU I I) ) ) 

4 SLMI I ) = VT 



MCE = 1 

number of constraint EVALUATIONS 

I = 1 

IF iKc I VII , 1) ) , EQ.O) GO TO 5 
I F (NPR .LE.O) GO TO 12 
WRITE (6,50) 

GO TO 12 

5 NFS = 1 

NUM3ER OF VERTICES IK) =2 "^TMES NO. CF VARIABLES. 

K = 2*NV 

number of DI SPLACEMENTS ALLOWED. 

NLIM = 5*NV+10 

NUMBER 0’= CONSECUTIVE TRIALS WITH UNCHANGED F£ TO 
T ERMIN AT E . 

NCT = NLIM+NV 
A L°HA=1 . 3D 0 
FK = K 
FKM=FK-1 .DO 
3ETA=ALPHA+ l.DO 

INSURE SEED OF RANDOM NUMBER GENE>^AT05 IS ODD. 

IQR = R^t'l.OT 

IF (MO DI IQR ,2 ) .£3 .0 ) I QR= I 3R<- 10 1 

SET UP INITIAL VERTICES 
FUN(1 ) = F6(V( 1, 1 ) ) 

Y*^N.= FUN(l) 

6 FI = l.DC 
FUNOLO = FUi\(l) 

DO 15 1=2, K 
FI = FI+l.DO 
LIMT = 0 

7 LIMT = LlMT+1 

C END CALCULATION IF FEASIBLE CENTROIC CANNOT EE FGUNC. 

IF ( L IMT .G£ .NLIM ) GO TO 11 
C 

DO 8 J=1 ,NV 
C 

C RANCOM number GENERATOR (RANDU) 

I CR = I GR*65539 

IF (IQR.LT.O) I3R = I OR + 21 474326474-1 
RCX = I QR 

RCX = RQX-.4656613D-9 

V(J,I) = BL ( J )4-R3X^( au( J )-3L( J ) ) 

I F( (V ( J , I ) -i-,5 00 ) .GT.2147433647 . )WRITE( 6, 100 ) 
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c 

c 



c 

c 



c 

c 

c 



c 

c 



c 

c 

c 



c 

c 

c 

c 



c 



IF (IP.EQ.l )V<J, I ) = IDINT(V( J, I ) + .5D0) 

8 CONTINUE 

OC 10 L=1,NLIM 
NC5 = NCE+1 

IF (KE( va, I ) ) .£Q ,0) GO TO 13 
00 9 J=1,NV 

VT = ( V( J , I) KEN{ J ) ) *.500 

IF( ( VT + .500 ). GT. 2147483647. )WR IT P (6» 100 ) 

IF (IP.EO.l) VT=I0IiNT (VT+.5D0) 

V (Jr n = VT 

9 CONTINUE 

10 CONTINUE 

11 IF (NPR.LE.O) GO TO 12 
WRITE (6r51) I 

CALL BOUT (NTrNPT rNFErNCErNV, NVTrVrl fFUN.CENr I) 

12 lER = -1 
GC TO 48 

15 DO 14 J=lrNV 

SUN( J) = SUM(J ) -t-V ( J , I J 

14 C£N( J) s SUM( J) /FI 

TRY TO ASSURE FEASIBLE CENTROID FOR STARTING. 

NCE = NCS+1 

IF ( KE( CEN ) .EQ. 0) GO TO 60 
SUM(J) = SUM(J) -V( J, I ) 

GO TO 7 
60 NFE = NFE+1 

FUN( I ) = FE (V(l , I n 

15 CONTINUE 

END CF LOOP SETTING OF INITIAL COMPLEX. 

I = ( NPR.LE. 0) GO TO 17 

CALL BOUT (NTrNPT,NFE,NCE,NVrNVT»V,K»FUNrCEN, 0) 

RIND THE WORST VERTEX, THE 'J’TH. 

J = 1 

00 16 1=2, K 

IF (FUN (J ) .G£ .FUN ( I) ) GC ^0 16 
J = I 

16 CCNTT NUE 

BASIC LOOP. ELIMINATE EACh WORST VERTEX IN TURN. 
MUST BECOME NO LONGER WORST, NOT MERELY I^PRCVED. 
N£X T-^O-WORST VERTEX, THE ’JN'“H ONE. 

17 JN = 1 

IF (J.EQ.l) JN = 2 

DC 18 I =1 ,K 

IF ( I .EQ.J ) GO TO 13 

IF (FUN ( JN ) .GE.FUNI I ) ) GO tq is 

JN = I 

13 CONTINUE 



LIMT=MUM8ER of moves DURING THIS TR I ^L TOW ARC TrE 
CENTROID DUE TO FUNCTION VALUE. 

Llf^ = 1 

COMPUTE centroid AMD OVER REFLECT WCRS'!' VER'EX. 

DO 19 I =l ,NV 
V^ = V( I , J) 

SUM( I ) = SU'^m-VT 

CEN( I ) = SUM (I) /PKM 

VT = 5ETA*CEN( I )-ALPHA = VT 

I F( (VT^-.500 ).GT .2147 433 647 . IWRITE (6,100 ) 

IF (IP.EQ.l) VT =I DINK VT+.5D0 ) 



IT 
F IN 
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c 

C INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 

19 V(I,J) = DMAXKOMINl (VT.eU (I) ) , BL <I ) ) 

C 

NT = NT+1 

CHECK FOR implicit CONSTRAINT VIOLATICN. 

20 DC 25 N=ltNLIM 
NCE = NCE+1 

IF (K E( V(1 , J ) ) .EQ .0) GO TO 26 

EVERY »KV*TH TI ME 1 0 VER-REFLECT THE CFFENOING VERTEX 
THROUGH THE BEST VERTEX. 

IF (MOO(N,KV) .NE.Ol GO TC 22 
CALL FBV (K, FUN, MI 

DC 21 I=1,NV 

VT = BETA*V( I tMI-ALPHA’^Vd ,J) ' 

IF( (VT + .5D0 ) .GT .2147 4326 47. ) HR I TE ( 6, 100) 

IF (IP.EQ.l) VT =101 NT IVT + .5D0 ) 

21 V(I,J) = DMAXK DMINKVT t 8U ( I ) ) ,8 L(I n 

GC TO 24 

CONSTRAINT VIOLATION: MOVE NEW POINT TOWARD CENTROID. 

22 DO 23 1=1, NV 
VT = (CEN( I ) +V (I, J ) I’^.SDO 
I = ( { VT^-, 5D0) .GT.214748 2647. ) WR IT = (6 ,100 ) 

IF (IP. £0.1) VT = IDIN T( VT^.SDO) 

V (I , J ) = VT 

23 CONTINUE 

2 A NT = MT+1 

25 CONTINUE 

I ER = 1 

CANNOT GET FEASIBLE VERTEX 8Y MOVING TOWARD CEN'RCID, 

OR 9Y OVER-REFLECTING THRU THE 8EST VERTEX. 

IF (MPR.LE. 0) GO TO 42 
WRITE ( 6, 52 ) NT, J 

CALL BOUT (NT,NPT ,NF£,MCE,NV,MVT,v,K,F’JN,C5M, J) 

GO TO 42. 

FEASIBLE vertex FOUND, EVALUATE THE OBJECTIVE !=UNCTiGN. 

26 NFE = NFE+1 
FUMTRY = FE( V( 1, J ) ) 

TEST TO SEE IF FUNCTION VALUE HAS NCT CHANGED. 

AFO =0ABS ( PJNTR Y-FUNOLD ) 

A fX=OMAXl ( DABS (EP’i'FUNGLC ) , £p ) 

activate the FOLLOWING TWO STATEMENTS FOR DIAGNOSTIC 
PURPOSES ONLY. 

WRITE ( 6, 9F) J, AFO, AM X , FUNTR Y , FUMOLC , FUN ( J ) ,FLN ( JN) , 
♦NTFS , N 

99 FCRMAT (IX, 12,6=15.7,215) 

IF ( AFO. GT .AMX) GO TO 27 
NTFS = NTF5+1 
IF (NTFS.L T.NCT) GO TO 28 
I ER = 0 

IF (NPR.LE.O) GO TO 42 
WRITE (6,53) K 

CALL BOUT (NT,NPT ,NFE,NCE,MV,NVT,v, K,FUN,CEN, C) 

GO TO 42 

27 NTFS = 0 

IS THE NEW VER'^EX NO LCNC-SR WORST? 

28 IF ( FUNTRY .LT .FUN ( JN ) ) GO TO 34 
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C TRIAL VERTEX IS STILL WORST; ADJUST TOWARD CENTROID. 

C EVERY ’KV’TH TIME, QVER-R £ = LEC't ThE Cr^ENOING VERTEX 

C THROUGH the BEST VERTEX. 

LIMT = LIYT+1 

IF (MOD (LIMT, KV) . NE.O ) GO TO 30 
CALL FBV (K,FUN,.M) 

C 

DO 29 I =1,NV 

VT = BETA*V(I,M)-AL?HA*V(I ,J) 

I F ( (VTt. 5 CO ) .GT .2i47A83 647 . )W R I TE ( 6 , 100 ) 

IF (IP.EQ.l) VT =IDINT( VT+.5D0 ) 

29 V(I,J) = DMAXK DM INK VT,BU(I ) ) ,BL(I ) ) 

C 

GO TO 32 

20 DO 31 I=1,NV 

VT = (CEN( n-H/( I, J ) )=*.500 
IF( (VT+.300 ) .GT .2147483647 . >WR IT E( 6, 100 ) 

1'= (IP.EQ.l) VT =IDI NT( VT+.5O0) 

V ( I,J ) = VT 

31 CONTINUE 

32 IF (L IMT .LT .NLI.M ) GO TO 33 

CANNOT make THE • J’ TH VERTEX NO LONGER WORST 3Y 
DISPLACING TOWARD THE CENTROID OR BY OVER-REFLEC TING 
THROUGH THE BEST VERTEX. 
lER = 2 

IF (NPR .LE. 0) GO TO 42 
WRITE (6,52) NT, J 

CALL BOUT ( NT,NRT,NFE ,NCE,NV,NVT ,V,K,FUN ,CEN, J) 

GO TO 42 

33 NT = nT+1 
GO TO 20 

SUCCESS; we HAVE A REPLACEMENT FOR VERTEX J. 

34 FUN(J) = FUNTRY 
FUNOLD = FUNTRY 
NPT = NPT+l 

EVERY IOO'TH PERMISSIBLE TR I AL , R ECOMOtjTE CENTROID 
summation to avoid CREEPING ERROR. 

IF (Mno( NPT, 100 ) .N£ .0 ) GO TO 37 

DO 3 6 I=1,NV 
SUM( I ) = 0.00 

DO 35 N=1,K 

25 5 JM( I ) = SUM( I )+•'/( I ,N ) 

CEN( I ) = SUM( I ) /PK 

36 CONTINUE 

LC = 0 
GO TO 39 

37 DO 3 8 ! =1 ,MV 

38 SUM(I ) = SUM( I)+V( I,J) 

LC = J 

39 IF (NOR ,lE .0 ) GO TO 40 
IF (MOD(MPT ,NPR) . NE.O) GO ^0 40 

CALL BOU"^ (NT ,NPT,NFE, NC£,NV,NVT ,V,K,FUN,CEN, LC ) 

HAS THE MAXIMUM NUMBER CF TRIALS BEEN REACHED WIThCU" 

convergence? if not, go to new TPIAL. 

40 IF ( NT. GE. iNTA) GO TO 41 

NEXT-TO-WORS"^ VERTEX NOW BECOMES WORST. 

J = JN 
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GC TO 17 

41 lER = 3 

IF (MPR .GT .0 ) WR ITE ( 6t 54) 

C 

C COLLECTOR POINT FOR ALL ENDINGS. 

C 1) CANNOT develop FEASI9LE VERTEX. 

C 2) CANNOT DEVELOP A NC-LCMGEP-WORST VERTEX. 

C 3) FUNCTION VALUE UNCHANGED FOR K TRIALS. 

C 4) , L iM IT ON TR IAL3 REACHED. 

C 5) CANNOT PINO FEASIBLE VERTEX AT START. 

42 CONTINUE 
C 

c find best vertex. 

CALL FBV (K,FUN,M) 

IF ( IER.GE.3 ) GO TO 44 
C 

C restart if this solution is significantly better THAN 

C THE PREVIOUStOR IF THIS IS THE FIRST TRY. 

IF ( NPR.LE.O) go TO 43 
WRITE (6t53) ( M,YMN,FUN( M) ) 

43 IF (FUN(M) .GE.YMN ) GO TC 47 

IF( DABS ( FUN( M)-YMN) . L£ .0.*<AX: (EP, EF^Y^n) ) GO TO 47 

C GIVE IT another TRY UNLESS LIMIT ON TRIALS REACHED. 

44 y'-N = FUN(M) 

FUN(l) = FJN(M) 

C 

DO 45 I =ltNV 
C 5N ( I ) = V ( I , M ) 

SUM (I ) = V ( I ,M) 

45 V( I ,1) = V( I ,M) 

C 

DC ^6 I=lfNVT 

46 XS(I) = V( I ,M) 

C 

IF (IER.LT.3) go to 6 

47 IF ( MPR.LE. 0) GO to 48 

CALL BOUT (NT,NpT ,NFE,NCE,NV,NVT,v,k,FUN»V( ,- l) 
WRITE (6.56) FUN(M) 

48 RETURN 
C 

49 FORMAT ( 'GIN CEX AND DIRECTION OF OUTLYING VARIA3LE't 
* ' AT START' ,I 5) 

50 format { '0 IMOl ICIT CONSTRAINT VIOLATED AT START.', 

*0 X, 'DEAD. END. ' ) 

51 FOPMAT( • OCANMOT FIND FE A SI 3LE ' , 14 , 

*'TH vertex or CBITROIC at START.') 

52 format ('CAT TRIAL ',14.' CANNOT FIND FEASIBLE VERTEX', 

WHICH IS NO LONGER WOR ST' , I 4, i 5 X , • RE STA R T FROM BEST', 
»• VERT EX .’ ) 

53 =ORMAT (40HOFUNCTION HAS BEEN ALMOST UNCHANGED FCR 15, 
•*7H TR I ALS) 

54 FCRMAT (27H0LIMIT CN TRIALS EXCEECED. ) 

55 FORMAT! • 08 E ST VERTEX IS NO . ' , 1 3 , ' CLD min WAS ', 
♦E15.7, ' NEW MIN IS ' , E15.7) 

56 =GRMAT t'OMIN OBJECTIVE FJNCtIOM IS ',E15.7) 

100 FORMAT ( 'O' , 'TRU.NCATI ON ERROR IN BOXOLX') 

END 



SUBROUTINE FBV (X,FUN,m) 
IMPLICIT RE AL»3 ( A-H, 0-Z) 

REAL *3 FUN( 50) 

M = 1 
C 

DC 1 I=2,K 

IF (FUN( M) .L£ .FUN ( I ) ) GO TO 1 
M = I 

1 CONTINUE 
RETURN 
END 



lER = 
lER = 
lER = 
lER = 
lER = 
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c 

c 



SUBROUTINE BOUT ( N T» N P T , NF E ,NC E ,N V , NVT , V , K ,FN ,C » I K) 
implicit realms (A-H,0-ZJ 
REAL*8 V( 50,50) ,FN(50) ,C(25) 

WRITE (6,4) NT,NPT,NFE,NC£ 

0 C 1 I =1,K 

WRITE (6,5) FN( I) , ( V( J ,I ) ,J=1 ,NV) 

IF (NVT.LE.NV) GO TO 1 
MVP = NV+1 

WRITE (6,6) ( V( J, I ) , J=NVP,.NVT) 

CCNTINUe 



IF ( IK .ME.O ) GO TO 2 
(C(I ) ,I =l 



WRITE (6,7) 
RETURN 

IF ( IK.GE.O ) 
WRITE (6,3) 
RETURN 
WRITE (6,9) 
RETURN 



,NV) 



GO TO 3 
(C( I) ,1=1 ,NV) 

IK, (C( I ) , 1 = 1, NV ) 



TR lALS _ = ’ , I = , 



’ , I" 

’ ,I 5,4X, 

f 



AX 



^ FCRMATCONO. TOTAL 
*'N0. FEASIBLE TRIALS = 

*'M0. FUNCTION EVALUATIONS = ' , I 5, AX, 
=*NQ, CONSTRAINT EVALUATICNS = ’,15/, 
«'0 function VALUE’, 6X, 

INDEPENDENT V ARI ASL ES / C £? ENOEN T OR 
*’CGNSTRAI NTS’ ) 



IMPLIC IT’ 



o 

7 

3 

9 



FORMAT 
FORMAT 
forma T 
FORM A 
FORMAT 
END 



{ IH ,Eia.7,2X, 7E1A. 7/{ 21X,7E1A.7 ) ) 

(21X,7EL A.7) 

( IGHOCENTRQIO IIX ,7E14 . 7 / ( 21 X , 7E 1 A . 7 ) ) 

{’0 BEST VERTEX' ,7X,7E14. 7/(21X,7ElA. 7) ) 
(’OCENTROID LESS VX • , 1 2 , 2X , 7E1 4. 7/ ( 2lX, 7 514 . 7 ) ) 



C 

C A /♦ 

C 

//GO. SYS IN 



CARD 

DO * 
•? 



GOES MERE 



1 0 



0.0 




0.005 


1 5. C 








13 


5G 


1 2 










1 00.0 




0.4 


10.0 




1.0 Q.O 


STEP 


14 


5C 


1 2 










16.0 




0.6 ' 


4 .0 




1.0 0,0 


STEP 


18 














1 3 1 


■2 


1100.0 


3.0 




7. 63636363 




2 3 1 


4 


-400,0 


3 .0 




A. 3 




3 3 2 


5 


-32.0 


4 ,3 




7. 0 




4 3 2 


6 


24. C 


4.3 




5. 73333333 




5 11 


i 


100 .0 










6 12 


Q 


16.0 










7 2 3 


7 


1. 0 


0.0 








3 2 4 


8 


1 .0 


0 .0 








9 2 5 


7 


1.0 


0 .0 








10 2 6 


s 


1.0 


0.0 








11 2 7 


9 


1.0 


7 .0 








12 2 31 0 


1 .0 


4.3 








13 2 91 


1 

A 


1. 0 


12.0 








14 2101 


-7 


1 .0 


2.0 








15 112 


9 


2.0 










16 line 


4. 0 










17 111 




-1.0 










13 112 


-7 


-1 .0 










1.0 




0. 0 


STEP 


1 






1 .0 




0.0 


STEP 


2 
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10 
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